Nektar++
NekLinSys.hpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File: NekLinSys.hpp
4 //
5 // For more information, please see: http://www.nektar.info
6 //
7 // The MIT License
8 //
9 // Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
10 // Department of Aeronautics, Imperial College London (UK), and Scientific
11 // Computing and Imaging Institute, University of Utah (USA).
12 //
13 // Permission is hereby granted, free of charge, to any person obtaining a
14 // copy of this software and associated documentation files (the "Software"),
15 // to deal in the Software without restriction, including without limitation
16 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
17 // and/or sell copies of the Software, and to permit persons to whom the
18 // Software is furnished to do so, subject to the following conditions:
19 //
20 // The above copyright notice and this permission notice shall be included
21 // in all copies or substantial portions of the Software.
22 //
23 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
24 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
25 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
26 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
27 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
28 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
29 // DEALINGS IN THE SOFTWARE.
30 //
31 // Description:
32 //
33 ///////////////////////////////////////////////////////////////////////////////
34 
35 #ifndef NEKTAR_LIB_UTILITIES_LINEAR_ALGEBRA_NEK_LINSYS_HPP
36 #define NEKTAR_LIB_UTILITIES_LINEAR_ALGEBRA_NEK_LINSYS_HPP
37 
38 #include <LibUtilities/BasicConst/NektarUnivTypeDefs.hpp>
39 #include <LibUtilities/LinearAlgebra/Lapack.hpp>
40 #include <LibUtilities/LinearAlgebra/NekMatrix.hpp>
41 #include <LibUtilities/LinearAlgebra/NekVector.hpp>
42 #include <LibUtilities/LinearAlgebra/MatrixType.h>
43 #include <LibUtilities/BasicUtils/ConsistentObjectAccess.hpp>
44 #include <LibUtilities/BasicUtils/RawType.hpp>
45 #include <LibUtilities/LinearAlgebra/NekLinSys.hpp>
46 #include <iostream>
47 
48 #include <type_traits>
49 
50 #ifdef max
51 #undef max
52 #endif
53 
54 #ifdef min
55 #undef min
56 #endif
57 
58 namespace Nektar
59 {
60  template<typename DataType>
61  struct IsSharedPointer : public std::false_type {};
62 
63  template<typename DataType>
64  struct IsSharedPointer<std::shared_ptr<DataType> > : public std::true_type {};
65 
66  // The solving of the linear system is located in this class instead of in the LinearSystem
67  // class because XCode gcc 4.2 didn't compile it correctly when it was moved to the
68  // LinearSystem class.
69  struct LinearSystemSolver
70  {
71  template<typename BVectorType, typename XVectorType>
72  static void Solve(const BVectorType& b, XVectorType& x, MatrixStorage m_matrixType,
73  const Array<OneD, const int>& m_ipivot, unsigned int n,
74  const Array<OneD, const double>& A,
75  char m_transposeFlag, unsigned int m_numberOfSubDiagonals,
76  unsigned int m_numberOfSuperDiagonals)
77  {
78  switch(m_matrixType)
79  {
80  case eFULL:
81  {
82  x = b;
83  int info = 0;
84  Lapack::Dgetrs('N',n,1,A.get(),n,(int *)m_ipivot.get(),x.GetRawPtr(),n,info);
85  if( info < 0 )
86  {
87  std::string message = "ERROR: The " + std::to_string(-info) + "th parameter had an illegal parameter for dgetrs";
88  ASSERTL0(false, message.c_str());
89  }
90 
91  }
92  break;
93  case eDIAGONAL:
94  for(unsigned int i = 0; i < A.num_elements(); ++i)
95  {
96  x[i] = b[i]*A[i];
97  }
98  break;
99  case eUPPER_TRIANGULAR:
100  {
101  x = b;
102  int info = 0;
103  Lapack::Dtptrs('U', m_transposeFlag, 'N', n, 1, A.get(), x.GetRawPtr(), n, info);
104 
105  if( info < 0 )
106  {
107  std::string message = "ERROR: The " + std::to_string(-info) + "th parameter had an illegal parameter for dtrtrs";
108  ASSERTL0(false, message.c_str());
109  }
110  else if( info > 0 )
111  {
112  std::string message = "ERROR: The " + std::to_string(-info) + "th diagonal element of A is 0 for dtrtrs";
113  ASSERTL0(false, message.c_str());
114  }
115  }
116  break;
117  case eLOWER_TRIANGULAR:
118  {
119  x = b;
120  int info = 0;
121  Lapack::Dtptrs('L', m_transposeFlag, 'N', n, 1, A.get(), x.GetRawPtr(), n, info);
122 
123  if( info < 0 )
124  {
125  std::string message = "ERROR: The " + std::to_string(-info) + "th parameter had an illegal parameter for dtrtrs";
126  ASSERTL0(false, message.c_str());
127  }
128  else if( info > 0 )
129  {
130  std::string message = "ERROR: The " + std::to_string(-info) + "th diagonal element of A is 0 for dtrtrs";
131  ASSERTL0(false, message.c_str());
132  }
133  }
134  break;
135  case eSYMMETRIC:
136  {
137  x = b;
138  int info = 0;
139  Lapack::Dsptrs('U', n, 1, A.get(), m_ipivot.get(), x.GetRawPtr(), x.GetRows(), info);
140  if( info < 0 )
141  {
142  std::string message = "ERROR: The " + std::to_string(-info) + "th parameter had an illegal parameter for dsptrs";
143  ASSERTL0(false, message.c_str());
144  }
145  }
146  break;
147  case ePOSITIVE_DEFINITE_SYMMETRIC:
148  {
149  x = b;
150  int info = 0;
151  Lapack::Dpptrs('U', n, 1, A.get(), x.GetRawPtr(), x.GetRows(), info);
152  if( info < 0 )
153  {
154  std::string message = "ERROR: The " + std::to_string(-info) + "th parameter had an illegal parameter for dpptrs";
155  ASSERTL0(false, message.c_str());
156  }
157  }
158  break;
159  case eBANDED:
160  {
161  x = b;
162  int KL = m_numberOfSubDiagonals;
163  int KU = m_numberOfSuperDiagonals;
164  int info = 0;
165 
166  Lapack::Dgbtrs(m_transposeFlag, n, KL, KU, 1, A.get(), 2*KL+KU+1, m_ipivot.get(), x.GetRawPtr(), n, info);
167 
168  if( info < 0 )
169  {
170  std::string message = "ERROR: The " + std::to_string(-info) + "th parameter had an illegal parameter for dgbtrs";
171  ASSERTL0(false, message.c_str());
172  }
173  }
174  break;
175  case ePOSITIVE_DEFINITE_SYMMETRIC_BANDED:
176  {
177  x = b;
178  int KU = m_numberOfSuperDiagonals;
179  int info = 0;
180 
181  Lapack::Dpbtrs('U', n, KU, 1, A.get(), KU+1, x.GetRawPtr(), n, info);
182 
183  if( info < 0 )
184  {
185  std::string message = "ERROR: The " + std::to_string(-info) + "th parameter had an illegal parameter for dpbtrs";
186  ASSERTL0(false, message.c_str());
187  }
188  }
189  break;
190  case eSYMMETRIC_BANDED:
191  NEKERROR(ErrorUtil::efatal, "Unhandled matrix type");
192  break;
193  case eUPPER_TRIANGULAR_BANDED:
194  NEKERROR(ErrorUtil::efatal, "Unhandled matrix type");
195  break;
196  case eLOWER_TRIANGULAR_BANDED:
197  NEKERROR(ErrorUtil::efatal, "Unhandled matrix type");
198  break;
199 
200  default:
201  NEKERROR(ErrorUtil::efatal, "Unhandled matrix type");
202  }
203 
204  }
205 
206  template<typename BVectorType, typename XVectorType>
207  static void SolveTranspose(const BVectorType& b, XVectorType& x, MatrixStorage m_matrixType,
208  const Array<OneD, const int>& m_ipivot, unsigned int n,
209  const Array<OneD, const double>& A,
210  char m_transposeFlag, unsigned int m_numberOfSubDiagonals,
211  unsigned int m_numberOfSuperDiagonals)
212  {
213  switch(m_matrixType)
214  {
215  case eFULL:
216  {
217  x = b;
218  int info = 0;
219  Lapack::Dgetrs('T',n,1,A.get(),n,(int *)m_ipivot.get(),x.GetRawPtr(), n,info);
220 
221  if( info < 0 )
222  {
223  std::string message = "ERROR: The " + std::to_string(-info) + "th parameter had an illegal parameter for dgetrs";
224  ASSERTL0(false, message.c_str());
225  }
226  }
227 
228  break;
229  case eDIAGONAL:
230  Solve(b, x, m_matrixType, m_ipivot, n, A, m_transposeFlag, m_numberOfSubDiagonals, m_numberOfSuperDiagonals);
231  break;
232  case eUPPER_TRIANGULAR:
233  {
234  char trans = m_transposeFlag;
235  if( trans == 'N' )
236  {
237  trans = 'T';
238  }
239  else
240  {
241  trans = 'N';
242  }
243 
244  x = b;
245  int info = 0;
246  Lapack::Dtptrs('U', trans, 'N', n, 1, A.get(), x.GetRawPtr(), n, info);
247 
248  if( info < 0 )
249  {
250  std::string message = "ERROR: The " + std::to_string(-info) + "th parameter had an illegal parameter for dtrtrs";
251  ASSERTL0(false, message.c_str());
252  }
253  else if( info > 0 )
254  {
255  std::string message = "ERROR: The " + std::to_string(-info) + "th diagonal element of A is 0 for dtrtrs";
256  ASSERTL0(false, message.c_str());
257  }
258  }
259 
260  break;
261  case eLOWER_TRIANGULAR:
262  {
263  char trans = m_transposeFlag;
264  if( trans == 'N' )
265  {
266  trans = 'T';
267  }
268  else
269  {
270  trans = 'N';
271  }
272  x = b;
273  int info = 0;
274  Lapack::Dtptrs('L', trans, 'N', n, 1, A.get(), x.GetRawPtr(), n, info);
275 
276  if( info < 0 )
277  {
278  std::string message = "ERROR: The " + std::to_string(-info) + "th parameter had an illegal parameter for dtrtrs";
279  ASSERTL0(false, message.c_str());
280  }
281  else if( info > 0 )
282  {
283  std::string message = "ERROR: The " + std::to_string(-info) + "th diagonal element of A is 0 for dtrtrs";
284  ASSERTL0(false, message.c_str());
285  }
286  }
287  break;
288  case eSYMMETRIC:
289  case ePOSITIVE_DEFINITE_SYMMETRIC:
290  case ePOSITIVE_DEFINITE_SYMMETRIC_BANDED:
291  Solve(b, x, m_matrixType, m_ipivot, n, A, m_transposeFlag, m_numberOfSubDiagonals, m_numberOfSuperDiagonals);
292  break;
293  case eBANDED:
294  {
295  x = b;
296  int KL = m_numberOfSubDiagonals;
297  int KU = m_numberOfSuperDiagonals;
298  int info = 0;
299 
300  Lapack::Dgbtrs(m_transposeFlag, n, KL, KU, 1, A.get(), 2*KL+KU+1, m_ipivot.get(), x.GetRawPtr(), n, info);
301 
302  if( info < 0 )
303  {
304  std::string message = "ERROR: The " + std::to_string(-info) + "th parameter had an illegal parameter for dgbtrs";
305  ASSERTL0(false, message.c_str());
306  }
307  }
308  break;
309  case eSYMMETRIC_BANDED:
310  NEKERROR(ErrorUtil::efatal, "Unhandled matrix type");
311  break;
312  case eUPPER_TRIANGULAR_BANDED:
313  NEKERROR(ErrorUtil::efatal, "Unhandled matrix type");
314  break;
315  case eLOWER_TRIANGULAR_BANDED:
316  NEKERROR(ErrorUtil::efatal, "Unhandled matrix type");
317  break;
318 
319  default:
320  NEKERROR(ErrorUtil::efatal, "Unhandled matrix type");
321  }
322  }
323  };
324 
325 
326  class LinearSystem
327  {
328  public:
329  template<typename MatrixType>
330  explicit LinearSystem(const std::shared_ptr<MatrixType> &theA, PointerWrapper wrapperType = eCopy) :
331  n(theA->GetRows()),
332  A(theA->GetPtr(), eVECTOR_WRAPPER),
333  m_ipivot(),
334  m_numberOfSubDiagonals(theA->GetNumberOfSubDiagonals()),
335  m_numberOfSuperDiagonals(theA->GetNumberOfSuperDiagonals()),
336  m_matrixType(theA->GetType()),
337  m_transposeFlag(theA->GetTransposeFlag())
338  {
339  // At some point we should fix this. We should upate the copy of
340  // A to be transposd for this to work.
341  ASSERTL0(theA->GetTransposeFlag() == 'N', "LinearSystem requires a non-transposed matrix.");
342  ASSERTL0( (wrapperType == eWrapper && theA->GetType() != eBANDED) || wrapperType == eCopy , "Banded matrices can't be wrapped");
343 
344  if( wrapperType == eCopy )
345  {
346  A = Array<OneD, double>(theA->GetPtr().num_elements());
347  CopyArray(theA->GetPtr(), A);
348  }
349 
350  FactorMatrix(*theA);
351  }
352 
353  template<typename MatrixType>
354  explicit LinearSystem(const MatrixType& theA, PointerWrapper wrapperType = eCopy) :
355  n(theA.GetRows()),
356  A(theA.GetPtr(), eVECTOR_WRAPPER),
357  m_ipivot(),
358  m_numberOfSubDiagonals(theA.GetNumberOfSubDiagonals()),
359  m_numberOfSuperDiagonals(theA.GetNumberOfSuperDiagonals()),
360  m_matrixType(theA.GetType()),
361  m_transposeFlag(theA.GetTransposeFlag())
362  {
363  // At some point we should fix this. We should upate the copy of
364  // A to be transposd for this to work.
365  ASSERTL0(theA.GetTransposeFlag() == 'N', "LinearSystem requires a non-transposed matrix.");
366  ASSERTL0( (wrapperType == eWrapper && theA.GetType() != eBANDED) || wrapperType == eCopy, "Banded matrices can't be wrapped" );
367 
368  if( wrapperType == eCopy )
369  {
370  A = Array<OneD, double>(theA.GetPtr().num_elements());
371  CopyArray(theA.GetPtr(), A);
372  }
373 
374  FactorMatrix(theA);
375  }
376 
377  LinearSystem(const LinearSystem& rhs) :
378  n(rhs.n),
379  A(rhs.A),
380  m_ipivot(rhs.m_ipivot),
381  m_numberOfSubDiagonals(rhs.m_numberOfSubDiagonals),
382  m_numberOfSuperDiagonals(rhs.m_numberOfSuperDiagonals),
383  m_matrixType(rhs.m_matrixType),
384  m_transposeFlag(rhs.m_transposeFlag)
385  {
386  }
387 
388  LinearSystem& operator=(const LinearSystem& rhs)
389  {
390  LinearSystem temp(rhs);
391  swap(temp);
392  return *this;
393  }
394 
395  ~LinearSystem() {}
396 
397  // In the following calls to Solve, VectorType must be a NekVector.
398  // Anything else won't compile.
399  template<typename VectorType>
400  RawType_t<VectorType> Solve(const VectorType& b)
401  {
402  RawType_t<VectorType> x(ConsistentObjectAccess<VectorType>::const_reference(b).GetRows());
403  LinearSystemSolver::Solve(ConsistentObjectAccess<VectorType>::const_reference(b), x, m_matrixType,
404  m_ipivot, n, A, m_transposeFlag, m_numberOfSubDiagonals, m_numberOfSuperDiagonals);
405  return x;
406  }
407 
408  template<typename BType, typename XType>
409  void Solve(const BType& b, XType& x) const
410  {
411  LinearSystemSolver::Solve(ConsistentObjectAccess<BType>::const_reference(b),
412  ConsistentObjectAccess<XType>::reference(x), m_matrixType,
413  m_ipivot, n, A, m_transposeFlag, m_numberOfSubDiagonals, m_numberOfSuperDiagonals);
414  }
415 
416  // Transpose variant of solve
417  template<typename VectorType>
418  RawType_t<VectorType> SolveTranspose(const VectorType& b)
419  {
420  RawType_t<VectorType> x(ConsistentObjectAccess<VectorType>::const_reference(b).GetRows());
421  LinearSystemSolver::SolveTranspose(ConsistentObjectAccess<VectorType>::const_reference(b), x, m_matrixType,
422  m_ipivot, n, A, m_transposeFlag, m_numberOfSubDiagonals, m_numberOfSuperDiagonals);
423  return x;
424  }
425 
426  template<typename BType, typename XType>
427  void SolveTranspose(const BType& b, XType& x) const
428  {
429  LinearSystemSolver::SolveTranspose(ConsistentObjectAccess<BType>::const_reference(b),
430  ConsistentObjectAccess<XType>::reference(x), m_matrixType,
431  m_ipivot, n, A, m_transposeFlag, m_numberOfSubDiagonals, m_numberOfSuperDiagonals);
432  }
433 
434  unsigned int GetRows() const { return n; }
435  unsigned int GetColumns() const { return n; }
436 
437  private:
438  template<typename MatrixType>
439  void FactorMatrix(const MatrixType& theA)
440  {
441  switch(m_matrixType)
442  {
443  case eFULL:
444  {
445  int m = theA.GetRows();
446  int n = theA.GetColumns();
447 
448  int pivotSize = std::max(1, std::min(m, n));
449  int info = 0;
450  m_ipivot = Array<OneD, int>(pivotSize);
451 
452  Lapack::Dgetrf(m, n, A.get(), m, m_ipivot.get(), info);
453 
454  if( info < 0 )
455  {
456  std::string message = "ERROR: The " + std::to_string(-info) + "th parameter had an illegal parameter for dgetrf";
457  ASSERTL0(false, message.c_str());
458  }
459  else if( info > 0 )
460  {
461  std::string message = "ERROR: Element u_" + std::to_string(info) + std::to_string(info) + " is 0 from dgetrf";
462  ASSERTL0(false, message.c_str());
463  }
464  }
465  break;
466  case eDIAGONAL:
467  for(unsigned int i = 0; i < theA.GetColumns(); ++i)
468  {
469  A[i] = 1.0/theA(i,i);
470  }
471  break;
472  case eUPPER_TRIANGULAR:
473  case eLOWER_TRIANGULAR:
474  break;
475  case eSYMMETRIC:
476  {
477  int info = 0;
478  int pivotSize = theA.GetRows();
479  m_ipivot = Array<OneD, int>(pivotSize);
480 
481  Lapack::Dsptrf('U', theA.GetRows(), A.get(), m_ipivot.get(), info);
482 
483  if( info < 0 )
484  {
485  std::string message = "ERROR: The " + std::to_string(-info) + "th parameter had an illegal parameter for dsptrf";
486  NEKERROR(ErrorUtil::efatal, message.c_str());
487  }
488  else if( info > 0 )
489  {
490  std::string message = "ERROR: Element u_" + std::to_string(info) + std::to_string(info) + " is 0 from dsptrf";
491  NEKERROR(ErrorUtil::efatal, message.c_str());
492  }
493  }
494  break;
495  case ePOSITIVE_DEFINITE_SYMMETRIC:
496  {
497  int info = 0;
498  Lapack::Dpptrf('U', theA.GetRows(), A.get(), info);
499 
500  if( info < 0 )
501  {
502  std::string message = "ERROR: The " + std::to_string(-info) + "th parameter had an illegal parameter for dpptrf";
503  NEKERROR(ErrorUtil::efatal, message.c_str());
504  }
505  else if( info > 0 )
506  {
507  std::string message = "ERROR: The leading minor of order " + std::to_string(info) + " is not positive definite from dpptrf";
508  NEKERROR(ErrorUtil::efatal, message.c_str());
509  }
510  }
511  break;
512  case eBANDED:
513  {
514  int M = n;
515  int N = n;
516  int KL = m_numberOfSubDiagonals;
517  int KU = m_numberOfSuperDiagonals;
518 
519  // The array we pass in to dgbtrf must have enough space for KL
520  // subdiagonals and KL+KU superdiagonals (see lapack users guide,
521  // in the section discussing band storage.
522  unsigned int requiredStorageSize = BandedMatrixFuncs::
523  GetRequiredStorageSize(n, n, KL, KL+KU);
524 
525  unsigned int rawRows = KL+KU+1;
526  A = Array<OneD, double>(requiredStorageSize);
527 
528  // Put the extra elements up front.
529  for(unsigned int i = 0; i < theA.GetColumns(); ++i)
530  {
531  std::copy(theA.GetRawPtr() + i*rawRows, theA.GetRawPtr() + (i+1)*rawRows,
532  A.get() + (i+1)*KL + i*rawRows);
533  }
534 
535  int info = 0;
536  int pivotSize = theA.GetRows();
537  m_ipivot = Array<OneD, int>(pivotSize);
538 
539  Lapack::Dgbtrf(M, N, KL, KU, A.get(), 2*KL+KU+1, m_ipivot.get(), info);
540 
541  if( info < 0 )
542  {
543  std::string message = "ERROR: The " + std::to_string(-info) + "th parameter had an illegal parameter for dgbtrf";
544  NEKERROR(ErrorUtil::efatal, message.c_str());
545  }
546  else if( info > 0 )
547  {
548  std::string message = "ERROR: Element u_" + std::to_string(info) + std::to_string(info) + " is 0 from dgbtrf";
549  NEKERROR(ErrorUtil::efatal, message.c_str());
550  }
551  }
552  break;
553  case ePOSITIVE_DEFINITE_SYMMETRIC_BANDED:
554  {
555  ASSERTL1(m_numberOfSuperDiagonals==m_numberOfSuperDiagonals,
556  std::string("Number of sub- and superdiagonals should ") +
557  std::string("be equal for a symmetric banded matrix"));
558 
559  int KU = m_numberOfSuperDiagonals;
560  int info = 0;
561  Lapack::Dpbtrf('U', theA.GetRows(), KU, A.get(), KU+1, info);
562 
563  if( info < 0 )
564  {
565  std::string message = "ERROR: The " + std::to_string(-info) + "th parameter had an illegal parameter for dpbtrf";
566  NEKERROR(ErrorUtil::efatal, message.c_str());
567  }
568  else if( info > 0 )
569  {
570  std::string message = "ERROR: The leading minor of order " + std::to_string(info) + " is not positive definite from dpbtrf";
571  NEKERROR(ErrorUtil::efatal, message.c_str());
572  }
573  }
574  break;
575  case eSYMMETRIC_BANDED:
576  NEKERROR(ErrorUtil::efatal, "Unhandled matrix type");
577  break;
578  case eUPPER_TRIANGULAR_BANDED:
579  NEKERROR(ErrorUtil::efatal, "Unhandled matrix type");
580  break;
581  case eLOWER_TRIANGULAR_BANDED:
582  NEKERROR(ErrorUtil::efatal, "Unhandled matrix type");
583  break;
584 
585  default:
586  NEKERROR(ErrorUtil::efatal, "Unhandled matrix type");
587  }
588  }
589 
590  void swap(LinearSystem& rhs)
591  {
592  std::swap(n, rhs.n);
593  std::swap(A, rhs.A);
594  std::swap(m_ipivot,rhs.m_ipivot);
595  std::swap(m_numberOfSubDiagonals, rhs.m_numberOfSubDiagonals);
596  std::swap(m_numberOfSuperDiagonals, rhs.m_numberOfSuperDiagonals);
597  std::swap(m_matrixType, rhs.m_matrixType);
598  std::swap(m_transposeFlag, rhs.m_transposeFlag);
599  }
600 
601  unsigned int n;
602  Array<OneD, double> A;
603  Array<OneD, int> m_ipivot;
604  unsigned int m_numberOfSubDiagonals;
605  unsigned int m_numberOfSuperDiagonals;
606  MatrixStorage m_matrixType;
607  char m_transposeFlag;
608  };
609 }
610 
611 #endif //NEKTAR_LIB_UTILITIES_LINEAR_ALGEBRA_NEK_LINSYS_HPP
Nektar::LinearSystem::SolveTranspose
void SolveTranspose(const BType &b, XType &x) const
Definition: NekLinSys.hpp:427
Lapack::Dgetrf
static void Dgetrf(const int &m, const int &n, double *a, const int &lda, int *ipiv, int &info)
General matrix LU factorisation.
Definition: Lapack.hpp:183
ASSERTL0
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:216
Nektar::LinearSystem::GetRows
unsigned int GetRows() const
Definition: NekLinSys.hpp:434
NEKERROR
#define NEKERROR(type, msg)
Assert Level 0 – Fundamental assert which is used whether in FULLDEBUG, DEBUG or OPT compilation mod...
Definition: ErrorUtil.hpp:209
Nektar::LinearSystem::A
Array< OneD, double > A
Definition: NekLinSys.hpp:602
Nektar::LinearSystem::m_matrixType
MatrixStorage m_matrixType
Definition: NekLinSys.hpp:606
Nektar::LinearSystem::Solve
RawType_t< VectorType > Solve(const VectorType &b)
Definition: NekLinSys.hpp:400
Nektar::LinearSystem::swap
void swap(LinearSystem &rhs)
Definition: NekLinSys.hpp:590
Nektar::LinearSystem::LinearSystem
LinearSystem(const std::shared_ptr< MatrixType > &theA, PointerWrapper wrapperType=eCopy)
Definition: NekLinSys.hpp:330
Lapack::Dsptrs
static void Dsptrs(const char &uplo, const int &n, const int &nrhs, const double *ap, const int *ipiv, double *b, const int &ldb, int &info)
Solve a real symmetric matrix problem using Bunch-Kaufman pivoting.
Definition: Lapack.hpp:115
std
STL namespace.
Nektar::LinearSystem::operator=
LinearSystem & operator=(const LinearSystem &rhs)
Definition: NekLinSys.hpp:388
Lapack::Dpbtrs
static void Dpbtrs(const char &uplo, const int &n, const int &kd, const int &nrhs, const double *ab, const int &ldab, double *b, const int &ldb, int &info)
Solve a real, Positive definite banded symmetric matrix problem using Cholesky factorization.
Definition: Lapack.hpp:157
CellMLToNektar.pycml.copy
def copy(self)
Definition: pycml.py:2663
NekVector.hpp
Lapack::Dgetrs
static void Dgetrs(const char &trans, const int &n, const int &nrhs, const double *a, const int &lda, int *ipiv, double *b, const int &ldb, int &info)
General matrix LU backsolve.
Definition: Lapack.hpp:190
RawType.hpp
Nektar::LinearSystem::m_numberOfSuperDiagonals
unsigned int m_numberOfSuperDiagonals
Definition: NekLinSys.hpp:605
Xxt::Solve
static void Solve(Nektar::Array< OneD, NekDouble > pX, struct crs_data *pCrs, Nektar::Array< OneD, NekDouble > pB)
Solve the matrix system for a given input vector b.
Definition: Xxt.hpp:188
Nektar::LinearSystem::GetColumns
unsigned int GetColumns() const
Definition: NekLinSys.hpp:435
Lapack::Dtptrs
static void Dtptrs(const char &uplo, const char &trans, const char &diag, const int &n, const int &nrhs, const double *a, double *b, const int &ldb, int &info)
Definition: Lapack.hpp:245
NektarUnivTypeDefs.hpp
Lapack::Dpbtrf
static void Dpbtrf(const char &uplo, const int &n, const int &kd, double *ab, const int &ldab, int &info)
Cholesky factorize a real positive-definite banded-symmetric matrix.
Definition: Lapack.hpp:148
Lapack.hpp
NekMatrix.hpp
Nektar::LinearSystem::SolveTranspose
RawType_t< VectorType > SolveTranspose(const VectorType &b)
Definition: NekLinSys.hpp:418
Nektar::LinearSystem::FactorMatrix
void FactorMatrix(const MatrixType &theA)
Definition: NekLinSys.hpp:439
Nektar::CopyArray
void CopyArray(const Array< OneD, ConstDataType > &source, Array< OneD, DataType > &dest)
Definition: SharedArray.hpp:765
Nektar::PointerWrapper
PointerWrapper
Specifies if the pointer passed to a NekMatrix or NekVector is copied into an internal representation...
Definition: PointerWrapper.h:42
Lapack::Dgbtrf
static void Dgbtrf(const int &m, const int &n, const int &kl, const int &ku, double *a, const int &lda, int *ipiv, int &info)
General banded matrix LU factorisation.
Definition: Lapack.hpp:166
Lapack::Dpptrf
static void Dpptrf(const char &uplo, const int &n, double *ap, int &info)
Cholesky factor a real Positive Definite packed-symmetric matrix.
Definition: Lapack.hpp:131
Nektar::LinearSystem::LinearSystem
LinearSystem(const LinearSystem &rhs)
Definition: NekLinSys.hpp:377
Lapack::Dpptrs
static void Dpptrs(const char &uplo, const int &n, const int &nrhs, const double *ap, double *b, const int &ldb, int &info)
Solve a real Positive defiinte symmetric matrix problem using Cholesky factorization.
Definition: Lapack.hpp:139
Lapack::Dsptrf
static void Dsptrf(const char &uplo, const int &n, double *ap, int *ipiv, int &info)
factor a real packed-symmetric matrix using Bunch-Kaufman pivoting.
Definition: Lapack.hpp:107
NekLinSys.hpp
MatrixType.h
ConsistentObjectAccess.hpp
Nektar::LinearSystem::LinearSystem
LinearSystem(const MatrixType &theA, PointerWrapper wrapperType=eCopy)
Definition: NekLinSys.hpp:354
Nektar::LinearSystemSolver::SolveTranspose
static void SolveTranspose(const BVectorType &b, XVectorType &x, MatrixStorage m_matrixType, const Array< OneD, const int > &m_ipivot, unsigned int n, const Array< OneD, const double > &A, char m_transposeFlag, unsigned int m_numberOfSubDiagonals, unsigned int m_numberOfSuperDiagonals)
Definition: NekLinSys.hpp:207
Nektar::LinearSystemSolver::Solve
static void Solve(const BVectorType &b, XVectorType &x, MatrixStorage m_matrixType, const Array< OneD, const int > &m_ipivot, unsigned int n, const Array< OneD, const double > &A, char m_transposeFlag, unsigned int m_numberOfSubDiagonals, unsigned int m_numberOfSuperDiagonals)
Definition: NekLinSys.hpp:72
Nektar::LinearSystem::Solve
void Solve(const BType &b, XType &x) const
Definition: NekLinSys.hpp:409
Nektar::rhs
StandardMatrixTag boost::call_traits< LhsDataType >::const_reference rhs
Definition: MatrixOperations.cpp:97
Nektar::LinearSystem::m_numberOfSubDiagonals
unsigned int m_numberOfSubDiagonals
Definition: NekLinSys.hpp:604
ASSERTL1
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode...
Definition: ErrorUtil.hpp:250
Lapack::Dgbtrs
static void Dgbtrs(const char &trans, const int &n, const int &kl, const int &ku, const int &nrhs, const double *a, const int &lda, const int *ipiv, double *b, const int &ldb, int &info)
Solve general banded matrix using LU factorisation.
Definition: Lapack.hpp:174
Nektar::RawType_t
typename RawType< T >::type RawType_t
Definition: RawType.hpp:76
Nektar::BandedMatrixFuncs::GetRequiredStorageSize
static unsigned int GetRequiredStorageSize(unsigned int totalRows, unsigned int totalColumns, unsigned int subDiags, unsigned int superDiags)
Calculates and returns the storage size required.
Definition: MatrixFuncs.cpp:44
Nektar::LinearSystem::m_ipivot
Array< OneD, int > m_ipivot
Definition: NekLinSys.hpp:603
Generated on Wed Jan 20 2021 03:33:37 for Nektar++ by   doxygen 1.8.13