Nektar++
Public Member Functions | Private Member Functions | Private Attributes | List of all members
Nektar::LinearSystem Class Reference

#include <NekLinSys.hpp>

Public Member Functions

template<typename MatrixType >
 LinearSystem (const std::shared_ptr< MatrixType > &theA, PointerWrapper wrapperType=eCopy)
 
template<typename MatrixType >
 LinearSystem (const MatrixType &theA, PointerWrapper wrapperType=eCopy)
 
 LinearSystem (const LinearSystem &rhs)
 
LinearSystem & operator= (const LinearSystem &rhs)
 
 ~LinearSystem ()
 
template<typename VectorType >
RawType_t< VectorType > Solve (const VectorType &b)
 
template<typename BType , typename XType >
void Solve (const BType &b, XType &x) const
 
template<typename VectorType >
RawType_t< VectorType > SolveTranspose (const VectorType &b)
 
template<typename BType , typename XType >
void SolveTranspose (const BType &b, XType &x) const
 
unsigned int GetRows () const
 
unsigned int GetColumns () const
 

Private Member Functions

template<typename MatrixType >
void FactorMatrix (const MatrixType &theA)
 
void swap (LinearSystem &rhs)
 

Private Attributes

unsigned int n
 
Array< OneD, double > A
 
Array< OneD, int > m_ipivot
 
unsigned int m_numberOfSubDiagonals
 
unsigned int m_numberOfSuperDiagonals
 
MatrixStorage m_matrixType
 
char m_transposeFlag
 

Detailed Description

Definition at line 326 of file NekLinSys.hpp.

Constructor & Destructor Documentation

◆ LinearSystem() [1/3]

template<typename MatrixType >
Nektar::LinearSystem::LinearSystem ( const std::shared_ptr< MatrixType > &  theA,
PointerWrapper  wrapperType = eCopy 
)
inlineexplicit

Definition at line 330 of file NekLinSys.hpp.

330  :
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().size());
347  CopyArray(theA->GetPtr(), A);
348  }
349 
350  FactorMatrix(*theA);
351  }
ASSERTL0
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:216
Nektar::LinearSystem::A
Array< OneD, double > A
Definition: NekLinSys.hpp:602
Nektar::LinearSystem::m_ipivot
Array< OneD, int > m_ipivot
Definition: NekLinSys.hpp:603
Nektar::LinearSystem::FactorMatrix
void FactorMatrix(const MatrixType &theA)
Definition: NekLinSys.hpp:439
Nektar::LinearSystem::m_numberOfSuperDiagonals
unsigned int m_numberOfSuperDiagonals
Definition: NekLinSys.hpp:605
Nektar::LinearSystem::m_numberOfSubDiagonals
unsigned int m_numberOfSubDiagonals
Definition: NekLinSys.hpp:604
Nektar::LinearSystem::m_matrixType
MatrixStorage m_matrixType
Definition: NekLinSys.hpp:606
Nektar::eVECTOR_WRAPPER
@ eVECTOR_WRAPPER
Definition: SharedArray.hpp:526
Nektar::CopyArray
void CopyArray(const Array< OneD, ConstDataType > &source, Array< OneD, DataType > &dest)
Definition: SharedArray.hpp:797

References ASSERTL0, Nektar::CopyArray(), Nektar::eBANDED, Nektar::eCopy, Nektar::eWrapper, and FactorMatrix().

◆ LinearSystem() [2/3]

template<typename MatrixType >
Nektar::LinearSystem::LinearSystem ( const MatrixType &  theA,
PointerWrapper  wrapperType = eCopy 
)
inlineexplicit

Definition at line 354 of file NekLinSys.hpp.

354  :
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().size());
371  CopyArray(theA.GetPtr(), A);
372  }
373 
374  FactorMatrix(theA);
375  }

References ASSERTL0, Nektar::CopyArray(), Nektar::eBANDED, Nektar::eCopy, Nektar::eWrapper, and FactorMatrix().

◆ LinearSystem() [3/3]

Nektar::LinearSystem::LinearSystem ( const LinearSystem &  rhs)
inline

Definition at line 377 of file NekLinSys.hpp.

377  :
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  }

◆ ~LinearSystem()

Nektar::LinearSystem::~LinearSystem ( )
inline

Definition at line 395 of file NekLinSys.hpp.

395 {}

Member Function Documentation

◆ FactorMatrix()

template<typename MatrixType >
void Nektar::LinearSystem::FactorMatrix ( const MatrixType &  theA)
inlineprivate

Definition at line 439 of file NekLinSys.hpp.

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  }
NEKERROR
#define NEKERROR(type, msg)
Assert Level 0 – Fundamental assert which is used whether in FULLDEBUG, DEBUG or OPT compilation mode...
Definition: ErrorUtil.hpp:209
ASSERTL1
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode....
Definition: ErrorUtil.hpp:250
CellMLToNektar.pycml.copy
def copy(self)
Definition: pycml.py:2663
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:274
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:305
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:239
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:182
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:256
Nektar::eLOWER_TRIANGULAR_BANDED
@ eLOWER_TRIANGULAR_BANDED
Definition: MatrixStorageType.h:53
Nektar::ePOSITIVE_DEFINITE_SYMMETRIC_BANDED
@ ePOSITIVE_DEFINITE_SYMMETRIC_BANDED
Definition: MatrixStorageType.h:51
Nektar::ePOSITIVE_DEFINITE_SYMMETRIC
@ ePOSITIVE_DEFINITE_SYMMETRIC
Definition: MatrixStorageType.h:48
Nektar::eUPPER_TRIANGULAR_BANDED
@ eUPPER_TRIANGULAR_BANDED
Definition: MatrixStorageType.h:52
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

References ASSERTL0, ASSERTL1, CellMLToNektar.pycml::copy(), Lapack::Dgbtrf(), Lapack::Dgetrf(), Lapack::Dpbtrf(), Lapack::Dpptrf(), Lapack::Dsptrf(), Nektar::eBANDED, Nektar::eDIAGONAL, Nektar::ErrorUtil::efatal, Nektar::eFULL, Nektar::eLOWER_TRIANGULAR, Nektar::eLOWER_TRIANGULAR_BANDED, Nektar::ePOSITIVE_DEFINITE_SYMMETRIC, Nektar::ePOSITIVE_DEFINITE_SYMMETRIC_BANDED, Nektar::eSYMMETRIC, Nektar::eSYMMETRIC_BANDED, Nektar::eUPPER_TRIANGULAR, Nektar::eUPPER_TRIANGULAR_BANDED, Nektar::BandedMatrixFuncs::GetRequiredStorageSize(), m_ipivot, m_matrixType, m_numberOfSubDiagonals, m_numberOfSuperDiagonals, n, and NEKERROR.

Referenced by LinearSystem().

◆ GetColumns()

unsigned int Nektar::LinearSystem::GetColumns ( ) const
inline

Definition at line 435 of file NekLinSys.hpp.

435 { return n; }

References n.

◆ GetRows()

unsigned int Nektar::LinearSystem::GetRows ( ) const
inline

Definition at line 434 of file NekLinSys.hpp.

434 { return n; }

References n.

Referenced by Solve(), and SolveTranspose().

◆ operator=()

LinearSystem& Nektar::LinearSystem::operator= ( const LinearSystem &  rhs)
inline

Definition at line 388 of file NekLinSys.hpp.

389  {
390  LinearSystem temp(rhs);
391  swap(temp);
392  return *this;
393  }
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

References swap().

◆ Solve() [1/2]

template<typename BType , typename XType >
void Nektar::LinearSystem::Solve ( const BType &  b,
XType &  x 
) const
inline

Definition at line 409 of file NekLinSys.hpp.

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  }
Nektar::ConsistentObjectAccess::const_reference
static const DataType & const_reference(const DataType &o)
Definition: ConsistentObjectAccess.hpp:50
Nektar::ConsistentObjectAccess::reference
static DataType & reference(DataType &o)
Definition: ConsistentObjectAccess.hpp:49
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

References m_ipivot, m_matrixType, m_numberOfSubDiagonals, m_numberOfSuperDiagonals, m_transposeFlag, n, and Nektar::LinearSystemSolver::Solve().

◆ Solve() [2/2]

template<typename VectorType >
RawType_t<VectorType> Nektar::LinearSystem::Solve ( const VectorType &  b)
inline

Definition at line 400 of file NekLinSys.hpp.

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  }
Nektar::LinearSystem::GetRows
unsigned int GetRows() const
Definition: NekLinSys.hpp:434

References GetRows(), m_ipivot, m_matrixType, m_numberOfSubDiagonals, m_numberOfSuperDiagonals, m_transposeFlag, n, and Nektar::LinearSystemSolver::Solve().

◆ SolveTranspose() [1/2]

template<typename BType , typename XType >
void Nektar::LinearSystem::SolveTranspose ( const BType &  b,
XType &  x 
) const
inline

Definition at line 427 of file NekLinSys.hpp.

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  }
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

References m_ipivot, m_matrixType, m_numberOfSubDiagonals, m_numberOfSuperDiagonals, m_transposeFlag, n, and Nektar::LinearSystemSolver::SolveTranspose().

◆ SolveTranspose() [2/2]

template<typename VectorType >
RawType_t<VectorType> Nektar::LinearSystem::SolveTranspose ( const VectorType &  b)
inline

Definition at line 418 of file NekLinSys.hpp.

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  }

References GetRows(), m_ipivot, m_matrixType, m_numberOfSubDiagonals, m_numberOfSuperDiagonals, m_transposeFlag, n, and Nektar::LinearSystemSolver::SolveTranspose().

Referenced by Nektar::LibUtilities::NodalUtil::GetWeights().

◆ swap()

void Nektar::LinearSystem::swap ( LinearSystem &  rhs)
inlineprivate

Definition at line 590 of file NekLinSys.hpp.

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  }

References A, m_ipivot, m_matrixType, m_numberOfSubDiagonals, m_numberOfSuperDiagonals, m_transposeFlag, and n.

Referenced by operator=().

Member Data Documentation

◆ A

Array<OneD, double> Nektar::LinearSystem::A
private

Definition at line 602 of file NekLinSys.hpp.

Referenced by swap().

◆ m_ipivot

Array<OneD, int> Nektar::LinearSystem::m_ipivot
private

Definition at line 603 of file NekLinSys.hpp.

Referenced by FactorMatrix(), Solve(), SolveTranspose(), and swap().

◆ m_matrixType

MatrixStorage Nektar::LinearSystem::m_matrixType
private

Definition at line 606 of file NekLinSys.hpp.

Referenced by FactorMatrix(), Solve(), SolveTranspose(), and swap().

◆ m_numberOfSubDiagonals

unsigned int Nektar::LinearSystem::m_numberOfSubDiagonals
private

Definition at line 604 of file NekLinSys.hpp.

Referenced by FactorMatrix(), Solve(), SolveTranspose(), and swap().

◆ m_numberOfSuperDiagonals

unsigned int Nektar::LinearSystem::m_numberOfSuperDiagonals
private

Definition at line 605 of file NekLinSys.hpp.

Referenced by FactorMatrix(), Solve(), SolveTranspose(), and swap().

◆ m_transposeFlag

char Nektar::LinearSystem::m_transposeFlag
private

Definition at line 607 of file NekLinSys.hpp.

Referenced by Solve(), SolveTranspose(), and swap().

◆ n

unsigned int Nektar::LinearSystem::n
private

Definition at line 601 of file NekLinSys.hpp.

Referenced by FactorMatrix(), GetColumns(), GetRows(), Solve(), SolveTranspose(), and swap().

Generated on Wed Nov 24 2021 10:33:30 for Nektar++ by doxygen 1.9.1