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>

Collaboration diagram for Nektar::LinearSystem:
Collaboration graph
[legend]

Public Member Functions

template<typename MatrixType >
 LinearSystem (const boost::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< VectorType >::type Solve (const VectorType &b)
 
template<typename BType , typename XType >
void Solve (const BType &b, XType &x) const
 
template<typename VectorType >
RawType< VectorType >::type 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 329 of file NekLinSys.hpp.

Constructor & Destructor Documentation

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

Definition at line 333 of file NekLinSys.hpp.

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

333  :
334  n(theA->GetRows()),
335  A(theA->GetPtr(), eVECTOR_WRAPPER),
336  m_ipivot(),
337  m_numberOfSubDiagonals(theA->GetNumberOfSubDiagonals()),
338  m_numberOfSuperDiagonals(theA->GetNumberOfSuperDiagonals()),
339  m_matrixType(theA->GetType()),
340  m_transposeFlag(theA->GetTransposeFlag())
341  {
342  // At some point we should fix this. We should upate the copy of
343  // A to be transposd for this to work.
344  ASSERTL0(theA->GetTransposeFlag() == 'N', "LinearSystem requires a non-transposed matrix.");
345  ASSERTL0( (wrapperType == eWrapper && theA->GetType() != eBANDED) || wrapperType == eCopy , "Banded matrices can't be wrapped");
346 
347  if( wrapperType == eCopy )
348  {
349  A = Array<OneD, double>(theA->GetPtr().num_elements());
350  CopyArray(theA->GetPtr(), A);
351  }
352 
353  FactorMatrix(*theA);
354  }
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:198
Array< OneD, double > A
Definition: NekLinSys.hpp:605
MatrixStorage m_matrixType
Definition: NekLinSys.hpp:609
unsigned int m_numberOfSuperDiagonals
Definition: NekLinSys.hpp:608
void FactorMatrix(const MatrixType &theA)
Definition: NekLinSys.hpp:442
void CopyArray(const Array< OneD, ConstDataType > &source, Array< OneD, DataType > &dest)
unsigned int m_numberOfSubDiagonals
Definition: NekLinSys.hpp:607
Array< OneD, int > m_ipivot
Definition: NekLinSys.hpp:606
template<typename MatrixType >
Nektar::LinearSystem::LinearSystem ( const MatrixType &  theA,
PointerWrapper  wrapperType = eCopy 
)
inlineexplicit

Definition at line 357 of file NekLinSys.hpp.

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

357  :
358  n(theA.GetRows()),
359  A(theA.GetPtr(), eVECTOR_WRAPPER),
360  m_ipivot(),
361  m_numberOfSubDiagonals(theA.GetNumberOfSubDiagonals()),
362  m_numberOfSuperDiagonals(theA.GetNumberOfSuperDiagonals()),
363  m_matrixType(theA.GetType()),
364  m_transposeFlag(theA.GetTransposeFlag())
365  {
366  // At some point we should fix this. We should upate the copy of
367  // A to be transposd for this to work.
368  ASSERTL0(theA.GetTransposeFlag() == 'N', "LinearSystem requires a non-transposed matrix.");
369  ASSERTL0( (wrapperType == eWrapper && theA.GetType() != eBANDED) || wrapperType == eCopy, "Banded matrices can't be wrapped" );
370 
371  if( wrapperType == eCopy )
372  {
373  A = Array<OneD, double>(theA.GetPtr().num_elements());
374  CopyArray(theA.GetPtr(), A);
375  }
376 
377  FactorMatrix(theA);
378  }
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:198
Array< OneD, double > A
Definition: NekLinSys.hpp:605
MatrixStorage m_matrixType
Definition: NekLinSys.hpp:609
unsigned int m_numberOfSuperDiagonals
Definition: NekLinSys.hpp:608
void FactorMatrix(const MatrixType &theA)
Definition: NekLinSys.hpp:442
void CopyArray(const Array< OneD, ConstDataType > &source, Array< OneD, DataType > &dest)
unsigned int m_numberOfSubDiagonals
Definition: NekLinSys.hpp:607
Array< OneD, int > m_ipivot
Definition: NekLinSys.hpp:606
Nektar::LinearSystem::LinearSystem ( const LinearSystem rhs)
inline

Definition at line 380 of file NekLinSys.hpp.

380  :
381  n(rhs.n),
382  A(rhs.A),
383  m_ipivot(rhs.m_ipivot),
384  m_numberOfSubDiagonals(rhs.m_numberOfSubDiagonals),
385  m_numberOfSuperDiagonals(rhs.m_numberOfSuperDiagonals),
386  m_matrixType(rhs.m_matrixType),
387  m_transposeFlag(rhs.m_transposeFlag)
388  {
389  }
Array< OneD, double > A
Definition: NekLinSys.hpp:605
MatrixStorage m_matrixType
Definition: NekLinSys.hpp:609
unsigned int m_numberOfSuperDiagonals
Definition: NekLinSys.hpp:608
unsigned int m_numberOfSubDiagonals
Definition: NekLinSys.hpp:607
Array< OneD, int > m_ipivot
Definition: NekLinSys.hpp:606
Nektar::LinearSystem::~LinearSystem ( )
inline

Definition at line 398 of file NekLinSys.hpp.

398 {}

Member Function Documentation

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

Definition at line 442 of file NekLinSys.hpp.

References ASSERTL0, ASSERTL1, CellMLToNektar.pycml::copy(), Nektar::eBANDED, Nektar::eDIAGONAL, 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().

443  {
444  switch(m_matrixType)
445  {
446  case eFULL:
447  {
448  int m = theA.GetRows();
449  int n = theA.GetColumns();
450 
451  int pivotSize = std::max(1, std::min(m, n));
452  int info = 0;
453  m_ipivot = Array<OneD, int>(pivotSize);
454 
455  Lapack::Dgetrf(m, n, A.get(), m, m_ipivot.get(), info);
456 
457  if( info < 0 )
458  {
459  std::string message = "ERROR: The " + boost::lexical_cast<std::string>(-info) + "th parameter had an illegal parameter for dgetrf";
460  ASSERTL0(false, message.c_str());
461  }
462  else if( info > 0 )
463  {
464  std::string message = "ERROR: Element u_" + boost::lexical_cast<std::string>(info) + boost::lexical_cast<std::string>(info) + " is 0 from dgetrf";
465  ASSERTL0(false, message.c_str());
466  }
467  }
468  break;
469  case eDIAGONAL:
470  for(unsigned int i = 0; i < theA.GetColumns(); ++i)
471  {
472  A[i] = 1.0/theA(i,i);
473  }
474  break;
475  case eUPPER_TRIANGULAR:
476  case eLOWER_TRIANGULAR:
477  break;
478  case eSYMMETRIC:
479  {
480  int info = 0;
481  int pivotSize = theA.GetRows();
482  m_ipivot = Array<OneD, int>(pivotSize);
483 
484  Lapack::Dsptrf('U', theA.GetRows(), A.get(), m_ipivot.get(), info);
485 
486  if( info < 0 )
487  {
488  std::string message = "ERROR: The " + boost::lexical_cast<std::string>(-info) + "th parameter had an illegal parameter for dsptrf";
489  ASSERTL0(false, message.c_str());
490  }
491  else if( info > 0 )
492  {
493  std::string message = "ERROR: Element u_" + boost::lexical_cast<std::string>(info) + boost::lexical_cast<std::string>(info) + " is 0 from dsptrf";
494  ASSERTL0(false, message.c_str());
495  }
496  }
497  break;
499  {
500  int info = 0;
501  Lapack::Dpptrf('U', theA.GetRows(), A.get(), info);
502 
503  if( info < 0 )
504  {
505  std::string message = "ERROR: The " + boost::lexical_cast<std::string>(-info) + "th parameter had an illegal parameter for dpptrf";
506  ASSERTL0(false, message.c_str());
507  }
508  else if( info > 0 )
509  {
510  std::string message = "ERROR: The leading minor of order " + boost::lexical_cast<std::string>(info) + " is not positive definite from dpptrf";
511  ASSERTL0(false, message.c_str());
512  }
513  }
514  break;
515  case eBANDED:
516  {
517  int M = n;
518  int N = n;
519  int KL = m_numberOfSubDiagonals;
520  int KU = m_numberOfSuperDiagonals;
521 
522  // The array we pass in to dgbtrf must have enough space for KL
523  // subdiagonals and KL+KU superdiagonals (see lapack users guide,
524  // in the section discussing band storage.
525  unsigned int requiredStorageSize = BandedMatrixFuncs::
526  GetRequiredStorageSize(n, n, KL, KL+KU);
527 
528  unsigned int rawRows = KL+KU+1;
529  A = Array<OneD, double>(requiredStorageSize);
530 
531  // Put the extra elements up front.
532  for(unsigned int i = 0; i < theA.GetColumns(); ++i)
533  {
534  std::copy(theA.GetRawPtr() + i*rawRows, theA.GetRawPtr() + (i+1)*rawRows,
535  A.get() + (i+1)*KL + i*rawRows);
536  }
537 
538  int info = 0;
539  int pivotSize = theA.GetRows();
540  m_ipivot = Array<OneD, int>(pivotSize);
541 
542  Lapack::Dgbtrf(M, N, KL, KU, A.get(), 2*KL+KU+1, m_ipivot.get(), info);
543 
544  if( info < 0 )
545  {
546  std::string message = "ERROR: The " + boost::lexical_cast<std::string>(-info) + "th parameter had an illegal parameter for dgbtrf";
547  ASSERTL0(false, message.c_str());
548  }
549  else if( info > 0 )
550  {
551  std::string message = "ERROR: Element u_" + boost::lexical_cast<std::string>(info) + boost::lexical_cast<std::string>(info) + " is 0 from dgbtrf";
552  ASSERTL0(false, message.c_str());
553  }
554  }
555  break;
557  {
559  std::string("Number of sub- and superdiagonals should ") +
560  std::string("be equal for a symmetric banded matrix"));
561 
562  int KU = m_numberOfSuperDiagonals;
563  int info = 0;
564  Lapack::Dpbtrf('U', theA.GetRows(), KU, A.get(), KU+1, info);
565 
566  if( info < 0 )
567  {
568  std::string message = "ERROR: The " + boost::lexical_cast<std::string>(-info) + "th parameter had an illegal parameter for dpbtrf";
569  ASSERTL0(false, message.c_str());
570  }
571  else if( info > 0 )
572  {
573  std::string message = "ERROR: The leading minor of order " + boost::lexical_cast<std::string>(info) + " is not positive definite from dpbtrf";
574  ASSERTL0(false, message.c_str());
575  }
576  }
577  break;
578  case eSYMMETRIC_BANDED:
579  NEKERROR(ErrorUtil::efatal, "Unhandled matrix type");
580  break;
582  NEKERROR(ErrorUtil::efatal, "Unhandled matrix type");
583  break;
585  NEKERROR(ErrorUtil::efatal, "Unhandled matrix type");
586  break;
587 
588  default:
589  NEKERROR(ErrorUtil::efatal, "Unhandled matrix type");
590  }
591  }
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:198
#define NEKERROR(type, msg)
Assert Level 0 – Fundamental assert which is used whether in FULLDEBUG, DEBUG or OPT compilation mod...
Definition: ErrorUtil.hpp:191
MatrixStorage m_matrixType
Definition: NekLinSys.hpp:609
unsigned int m_numberOfSuperDiagonals
Definition: NekLinSys.hpp:608
unsigned int m_numberOfSubDiagonals
Definition: NekLinSys.hpp:607
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode...
Definition: ErrorUtil.hpp:228
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:43
Array< OneD, int > m_ipivot
Definition: NekLinSys.hpp:606
unsigned int Nektar::LinearSystem::GetColumns ( ) const
inline

Definition at line 438 of file NekLinSys.hpp.

References n.

438 { return n; }
unsigned int Nektar::LinearSystem::GetRows ( ) const
inline

Definition at line 437 of file NekLinSys.hpp.

References n.

Referenced by Solve(), and SolveTranspose().

437 { return n; }
LinearSystem& Nektar::LinearSystem::operator= ( const LinearSystem rhs)
inline

Definition at line 391 of file NekLinSys.hpp.

References swap().

392  {
393  LinearSystem temp(rhs);
394  swap(temp);
395  return *this;
396  }
void swap(LinearSystem &rhs)
Definition: NekLinSys.hpp:593
LinearSystem(const boost::shared_ptr< MatrixType > &theA, PointerWrapper wrapperType=eCopy)
Definition: NekLinSys.hpp:333
template<typename VectorType >
RawType<VectorType>::type Nektar::LinearSystem::Solve ( const VectorType &  b)
inline

Definition at line 403 of file NekLinSys.hpp.

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

Referenced by Nektar::StdRegions::Invert().

404  {
405  typename RawType<VectorType>::type x(ConsistentObjectAccess<VectorType>::const_reference(b).GetRows());
406  LinearSystemSolver::Solve(ConsistentObjectAccess<VectorType>::const_reference(b), x, m_matrixType,
408  return x;
409  }
MatrixStorage m_matrixType
Definition: NekLinSys.hpp:609
unsigned int m_numberOfSuperDiagonals
Definition: NekLinSys.hpp:608
unsigned int GetRows() const
Definition: NekLinSys.hpp:437
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
unsigned int m_numberOfSubDiagonals
Definition: NekLinSys.hpp:607
Array< OneD, int > m_ipivot
Definition: NekLinSys.hpp:606
template<typename BType , typename XType >
void Nektar::LinearSystem::Solve ( const BType &  b,
XType &  x 
) const
inline

Definition at line 412 of file NekLinSys.hpp.

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

413  {
414  LinearSystemSolver::Solve(ConsistentObjectAccess<BType>::const_reference(b),
415  ConsistentObjectAccess<XType>::reference(x), m_matrixType,
417  }
MatrixStorage m_matrixType
Definition: NekLinSys.hpp:609
unsigned int m_numberOfSuperDiagonals
Definition: NekLinSys.hpp:608
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
unsigned int m_numberOfSubDiagonals
Definition: NekLinSys.hpp:607
Array< OneD, int > m_ipivot
Definition: NekLinSys.hpp:606
template<typename VectorType >
RawType<VectorType>::type Nektar::LinearSystem::SolveTranspose ( const VectorType &  b)
inline

Definition at line 421 of file NekLinSys.hpp.

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

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

422  {
423  typename RawType<VectorType>::type x(ConsistentObjectAccess<VectorType>::const_reference(b).GetRows());
424  LinearSystemSolver::SolveTranspose(ConsistentObjectAccess<VectorType>::const_reference(b), x, m_matrixType,
426  return x;
427  }
MatrixStorage m_matrixType
Definition: NekLinSys.hpp:609
unsigned int m_numberOfSuperDiagonals
Definition: NekLinSys.hpp:608
unsigned int GetRows() const
Definition: NekLinSys.hpp:437
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:210
unsigned int m_numberOfSubDiagonals
Definition: NekLinSys.hpp:607
Array< OneD, int > m_ipivot
Definition: NekLinSys.hpp:606
template<typename BType , typename XType >
void Nektar::LinearSystem::SolveTranspose ( const BType &  b,
XType &  x 
) const
inline

Definition at line 430 of file NekLinSys.hpp.

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

431  {
432  LinearSystemSolver::SolveTranspose(ConsistentObjectAccess<BType>::const_reference(b),
433  ConsistentObjectAccess<XType>::reference(x), m_matrixType,
435  }
MatrixStorage m_matrixType
Definition: NekLinSys.hpp:609
unsigned int m_numberOfSuperDiagonals
Definition: NekLinSys.hpp:608
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:210
unsigned int m_numberOfSubDiagonals
Definition: NekLinSys.hpp:607
Array< OneD, int > m_ipivot
Definition: NekLinSys.hpp:606
void Nektar::LinearSystem::swap ( LinearSystem rhs)
inlineprivate

Definition at line 593 of file NekLinSys.hpp.

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

Referenced by operator=().

594  {
595  std::swap(n, rhs.n);
596  std::swap(A, rhs.A);
597  std::swap(m_ipivot,rhs.m_ipivot);
598  std::swap(m_numberOfSubDiagonals, rhs.m_numberOfSubDiagonals);
599  std::swap(m_numberOfSuperDiagonals, rhs.m_numberOfSuperDiagonals);
600  std::swap(m_matrixType, rhs.m_matrixType);
601  std::swap(m_transposeFlag, rhs.m_transposeFlag);
602  }
MatrixStorage m_matrixType
Definition: NekLinSys.hpp:609
unsigned int m_numberOfSuperDiagonals
Definition: NekLinSys.hpp:608
unsigned int m_numberOfSubDiagonals
Definition: NekLinSys.hpp:607
Array< OneD, int > m_ipivot
Definition: NekLinSys.hpp:606

Member Data Documentation

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

Definition at line 605 of file NekLinSys.hpp.

Referenced by LinearSystem(), and swap().

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

Definition at line 606 of file NekLinSys.hpp.

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

MatrixStorage Nektar::LinearSystem::m_matrixType
private

Definition at line 609 of file NekLinSys.hpp.

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

unsigned int Nektar::LinearSystem::m_numberOfSubDiagonals
private

Definition at line 607 of file NekLinSys.hpp.

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

unsigned int Nektar::LinearSystem::m_numberOfSuperDiagonals
private

Definition at line 608 of file NekLinSys.hpp.

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

char Nektar::LinearSystem::m_transposeFlag
private

Definition at line 610 of file NekLinSys.hpp.

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

unsigned int Nektar::LinearSystem::n
private

Definition at line 604 of file NekLinSys.hpp.

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