Nektar++
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
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)
 
LinearSystemoperator= (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 375 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 379 of file NekLinSys.hpp.

381  : n(theA->GetRows()), A(theA->GetPtr(), eVECTOR_WRAPPER), m_ipivot(),
382  m_numberOfSubDiagonals(theA->GetNumberOfSubDiagonals()),
383  m_numberOfSuperDiagonals(theA->GetNumberOfSuperDiagonals()),
384  m_matrixType(theA->GetType()),
385  m_transposeFlag(theA->GetTransposeFlag())
386  {
387  // At some point we should fix this. We should upate the copy of
388  // A to be transposd for this to work.
389  ASSERTL0(theA->GetTransposeFlag() == 'N',
390  "LinearSystem requires a non-transposed matrix.");
391  ASSERTL0((wrapperType == eWrapper && theA->GetType() != eBANDED) ||
392  wrapperType == eCopy,
393  "Banded matrices can't be wrapped");
394 
395  if (wrapperType == eCopy)
396  {
397  A = Array<OneD, double>(theA->GetPtr().size());
398  CopyArray(theA->GetPtr(), A);
399  }
400 
401  FactorMatrix(*theA);
402  }
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:215
Array< OneD, double > A
Definition: NekLinSys.hpp:695
Array< OneD, int > m_ipivot
Definition: NekLinSys.hpp:696
void FactorMatrix(const MatrixType &theA)
Definition: NekLinSys.hpp:505
unsigned int m_numberOfSuperDiagonals
Definition: NekLinSys.hpp:698
unsigned int m_numberOfSubDiagonals
Definition: NekLinSys.hpp:697
MatrixStorage m_matrixType
Definition: NekLinSys.hpp:699
@ eVECTOR_WRAPPER
void CopyArray(const Array< OneD, ConstDataType > &source, Array< OneD, DataType > &dest)

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 405 of file NekLinSys.hpp.

407  : n(theA.GetRows()), A(theA.GetPtr(), eVECTOR_WRAPPER), m_ipivot(),
408  m_numberOfSubDiagonals(theA.GetNumberOfSubDiagonals()),
409  m_numberOfSuperDiagonals(theA.GetNumberOfSuperDiagonals()),
410  m_matrixType(theA.GetType()), m_transposeFlag(theA.GetTransposeFlag())
411  {
412  // At some point we should fix this. We should upate the copy of
413  // A to be transposd for this to work.
414  ASSERTL0(theA.GetTransposeFlag() == 'N',
415  "LinearSystem requires a non-transposed matrix.");
416  ASSERTL0((wrapperType == eWrapper && theA.GetType() != eBANDED) ||
417  wrapperType == eCopy,
418  "Banded matrices can't be wrapped");
419 
420  if (wrapperType == eCopy)
421  {
422  A = Array<OneD, double>(theA.GetPtr().size());
423  CopyArray(theA.GetPtr(), A);
424  }
425 
426  FactorMatrix(theA);
427  }

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 429 of file NekLinSys.hpp.

430  : n(rhs.n), A(rhs.A), m_ipivot(rhs.m_ipivot),
431  m_numberOfSubDiagonals(rhs.m_numberOfSubDiagonals),
432  m_numberOfSuperDiagonals(rhs.m_numberOfSuperDiagonals),
433  m_matrixType(rhs.m_matrixType), m_transposeFlag(rhs.m_transposeFlag)
434  {
435  }

◆ ~LinearSystem()

Nektar::LinearSystem::~LinearSystem ( )
inline

Definition at line 444 of file NekLinSys.hpp.

445  {
446  }

Member Function Documentation

◆ FactorMatrix()

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

Definition at line 505 of file NekLinSys.hpp.

506  {
507  switch (m_matrixType)
508  {
509  case eFULL:
510  {
511  int m = theA.GetRows();
512  int n = theA.GetColumns();
513 
514  int pivotSize = std::max(1, std::min(m, n));
515  int info = 0;
516  m_ipivot = Array<OneD, int>(pivotSize);
517 
518  Lapack::Dgetrf(m, n, A.get(), m, m_ipivot.get(), info);
519 
520  if (info < 0)
521  {
522  std::string message =
523  "ERROR: The " + std::to_string(-info) +
524  "th parameter had an illegal parameter for dgetrf";
525  ASSERTL0(false, message.c_str());
526  }
527  else if (info > 0)
528  {
529  std::string message =
530  "ERROR: Element u_" + std::to_string(info) +
531  std::to_string(info) + " is 0 from dgetrf";
532  ASSERTL0(false, message.c_str());
533  }
534  }
535  break;
536  case eDIAGONAL:
537  for (unsigned int i = 0; i < theA.GetColumns(); ++i)
538  {
539  A[i] = 1.0 / theA(i, i);
540  }
541  break;
542  case eUPPER_TRIANGULAR:
543  case eLOWER_TRIANGULAR:
544  break;
545  case eSYMMETRIC:
546  {
547  int info = 0;
548  int pivotSize = theA.GetRows();
549  m_ipivot = Array<OneD, int>(pivotSize);
550 
551  Lapack::Dsptrf('U', theA.GetRows(), A.get(), m_ipivot.get(),
552  info);
553 
554  if (info < 0)
555  {
556  std::string message =
557  "ERROR: The " + std::to_string(-info) +
558  "th parameter had an illegal parameter for dsptrf";
559  NEKERROR(ErrorUtil::efatal, message.c_str());
560  }
561  else if (info > 0)
562  {
563  std::string message =
564  "ERROR: Element u_" + std::to_string(info) +
565  std::to_string(info) + " is 0 from dsptrf";
566  NEKERROR(ErrorUtil::efatal, message.c_str());
567  }
568  }
569  break;
571  {
572  int info = 0;
573  Lapack::Dpptrf('U', theA.GetRows(), A.get(), info);
574 
575  if (info < 0)
576  {
577  std::string message =
578  "ERROR: The " + std::to_string(-info) +
579  "th parameter had an illegal parameter for dpptrf";
580  NEKERROR(ErrorUtil::efatal, message.c_str());
581  }
582  else if (info > 0)
583  {
584  std::string message =
585  "ERROR: The leading minor of order " +
586  std::to_string(info) +
587  " is not positive definite from dpptrf";
588  NEKERROR(ErrorUtil::efatal, message.c_str());
589  }
590  }
591  break;
592  case eBANDED:
593  {
594  int M = n;
595  int N = n;
596  int KL = m_numberOfSubDiagonals;
597  int KU = m_numberOfSuperDiagonals;
598 
599  // The array we pass in to dgbtrf must have enough space for KL
600  // subdiagonals and KL+KU superdiagonals (see lapack users
601  // guide, in the section discussing band storage.
602  unsigned int requiredStorageSize =
604  KL + KU);
605 
606  unsigned int rawRows = KL + KU + 1;
607  A = Array<OneD, double>(requiredStorageSize);
608 
609  // Put the extra elements up front.
610  for (unsigned int i = 0; i < theA.GetColumns(); ++i)
611  {
612  std::copy(theA.GetRawPtr() + i * rawRows,
613  theA.GetRawPtr() + (i + 1) * rawRows,
614  A.get() + (i + 1) * KL + i * rawRows);
615  }
616 
617  int info = 0;
618  int pivotSize = theA.GetRows();
619  m_ipivot = Array<OneD, int>(pivotSize);
620 
621  Lapack::Dgbtrf(M, N, KL, KU, A.get(), 2 * KL + KU + 1,
622  m_ipivot.get(), info);
623 
624  if (info < 0)
625  {
626  std::string message =
627  "ERROR: The " + std::to_string(-info) +
628  "th parameter had an illegal parameter for dgbtrf";
629  NEKERROR(ErrorUtil::efatal, message.c_str());
630  }
631  else if (info > 0)
632  {
633  std::string message =
634  "ERROR: Element u_" + std::to_string(info) +
635  std::to_string(info) + " is 0 from dgbtrf";
636  NEKERROR(ErrorUtil::efatal, message.c_str());
637  }
638  }
639  break;
641  {
642  ASSERTL1(
644  std::string("Number of sub- and superdiagonals should ") +
645  std::string("be equal for a symmetric banded matrix"));
646 
647  int KU = m_numberOfSuperDiagonals;
648  int info = 0;
649  Lapack::Dpbtrf('U', theA.GetRows(), KU, A.get(), KU + 1, info);
650 
651  if (info < 0)
652  {
653  std::string message =
654  "ERROR: The " + std::to_string(-info) +
655  "th parameter had an illegal parameter for dpbtrf";
656  NEKERROR(ErrorUtil::efatal, message.c_str());
657  }
658  else if (info > 0)
659  {
660  std::string message =
661  "ERROR: The leading minor of order " +
662  std::to_string(info) +
663  " is not positive definite from dpbtrf";
664  NEKERROR(ErrorUtil::efatal, message.c_str());
665  }
666  }
667  break;
668  case eSYMMETRIC_BANDED:
669  NEKERROR(ErrorUtil::efatal, "Unhandled matrix type");
670  break;
672  NEKERROR(ErrorUtil::efatal, "Unhandled matrix type");
673  break;
675  NEKERROR(ErrorUtil::efatal, "Unhandled matrix type");
676  break;
677 
678  default:
679  NEKERROR(ErrorUtil::efatal, "Unhandled matrix type");
680  }
681  }
#define NEKERROR(type, msg)
Assert Level 0 – Fundamental assert which is used whether in FULLDEBUG, DEBUG or OPT compilation mode...
Definition: ErrorUtil.hpp:209
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode....
Definition: ErrorUtil.hpp:249
def copy(self)
Definition: pycml.py:2663
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:250
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:281
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:218
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:165
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:234
@ eLOWER_TRIANGULAR_BANDED
@ ePOSITIVE_DEFINITE_SYMMETRIC_BANDED
@ ePOSITIVE_DEFINITE_SYMMETRIC
@ eUPPER_TRIANGULAR_BANDED
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 499 of file NekLinSys.hpp.

500  {
501  return n;
502  }

References n.

◆ GetRows()

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

Definition at line 495 of file NekLinSys.hpp.

496  {
497  return n;
498  }

References n.

Referenced by Solve(), and SolveTranspose().

◆ operator=()

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

Definition at line 437 of file NekLinSys.hpp.

438  {
439  LinearSystem temp(rhs);
440  swap(temp);
441  return *this;
442  }
void swap(LinearSystem &rhs)
Definition: NekLinSys.hpp:683
LinearSystem(const std::shared_ptr< MatrixType > &theA, PointerWrapper wrapperType=eCopy)
Definition: NekLinSys.hpp:379

References swap().

◆ Solve() [1/2]

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

Definition at line 463 of file NekLinSys.hpp.

464  {
470  }
static const DataType & const_reference(const DataType &o)
static DataType & reference(DataType &o)
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:75

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

◆ SolveTranspose() [1/2]

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

Definition at line 486 of file NekLinSys.hpp.

487  {
493  }
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:235

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

◆ swap()

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

Definition at line 683 of file NekLinSys.hpp.

684  {
685  std::swap(n, rhs.n);
686  std::swap(A, rhs.A);
687  std::swap(m_ipivot, rhs.m_ipivot);
688  std::swap(m_numberOfSubDiagonals, rhs.m_numberOfSubDiagonals);
689  std::swap(m_numberOfSuperDiagonals, rhs.m_numberOfSuperDiagonals);
690  std::swap(m_matrixType, rhs.m_matrixType);
691  std::swap(m_transposeFlag, rhs.m_transposeFlag);
692  }

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 695 of file NekLinSys.hpp.

Referenced by swap().

◆ m_ipivot

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

Definition at line 696 of file NekLinSys.hpp.

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

◆ m_matrixType

MatrixStorage Nektar::LinearSystem::m_matrixType
private

Definition at line 699 of file NekLinSys.hpp.

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

◆ m_numberOfSubDiagonals

unsigned int Nektar::LinearSystem::m_numberOfSubDiagonals
private

Definition at line 697 of file NekLinSys.hpp.

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

◆ m_numberOfSuperDiagonals

unsigned int Nektar::LinearSystem::m_numberOfSuperDiagonals
private

Definition at line 698 of file NekLinSys.hpp.

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

◆ m_transposeFlag

char Nektar::LinearSystem::m_transposeFlag
private

Definition at line 700 of file NekLinSys.hpp.

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

◆ n

unsigned int Nektar::LinearSystem::n
private

Definition at line 694 of file NekLinSys.hpp.

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