#include <NekLinSys.hpp>

Collaboration diagram for Nektar::LinearSystem:

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)

LinearSystem &	operator= (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().

                                                                                                                :
                 n(theA->GetRows()),
                 A(theA->GetPtr(), eVECTOR_WRAPPER),
                 m_ipivot(),
                 m_numberOfSubDiagonals(theA->GetNumberOfSubDiagonals()),
                 m_numberOfSuperDiagonals(theA->GetNumberOfSuperDiagonals()),
                 m_matrixType(theA->GetType()),
                 m_transposeFlag(theA->GetTransposeFlag())
             {
                 // At some point we should fix this.  We should upate the copy of 
                 // A to be transposd for this to work.
                 ASSERTL0(theA->GetTransposeFlag() == 'N', "LinearSystem requires a non-transposed matrix.");
                 ASSERTL0( (wrapperType == eWrapper && theA->GetType() != eBANDED) || wrapperType == eCopy , "Banded matrices can't be wrapped");
                 
                 if( wrapperType == eCopy )
                 {
                     A = Array<OneD, double>(theA->GetPtr().num_elements());
                     CopyArray(theA->GetPtr(), A);
                 }
                 
                 FactorMatrix(*theA);
             }

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().

                                                                                               :
                 n(theA.GetRows()),
                 A(theA.GetPtr(), eVECTOR_WRAPPER),
                 m_ipivot(),
                 m_numberOfSubDiagonals(theA.GetNumberOfSubDiagonals()),
                 m_numberOfSuperDiagonals(theA.GetNumberOfSuperDiagonals()),
                 m_matrixType(theA.GetType()),
                 m_transposeFlag(theA.GetTransposeFlag())
             {
                 // At some point we should fix this.  We should upate the copy of 
                 // A to be transposd for this to work.
                 ASSERTL0(theA.GetTransposeFlag() == 'N', "LinearSystem requires a non-transposed matrix.");
                 ASSERTL0( (wrapperType == eWrapper && theA.GetType() != eBANDED) || wrapperType == eCopy, "Banded matrices can't be wrapped" );
                 
                 if( wrapperType == eCopy )
                 {
                     A = Array<OneD, double>(theA.GetPtr().num_elements());
                     CopyArray(theA.GetPtr(), A);
                 }
 
                 FactorMatrix(theA);
             }

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

inline

Definition at line 380 of file NekLinSys.hpp.

                                                   :
                 n(rhs.n),
                 A(rhs.A),
                 m_ipivot(rhs.m_ipivot),
                 m_numberOfSubDiagonals(rhs.m_numberOfSubDiagonals),
                 m_numberOfSuperDiagonals(rhs.m_numberOfSuperDiagonals),
                 m_matrixType(rhs.m_matrixType),
                 m_transposeFlag(rhs.m_transposeFlag)
             {
             }

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 A, ASSERTL0, ASSERTL1, 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().

             {
                 switch(m_matrixType)
                 {
                     case eFULL:
                         {
                             int m = theA.GetRows();
                             int n = theA.GetColumns();
                             
                             int pivotSize = std::max(1, std::min(m, n));
                             int info = 0;
                             m_ipivot = Array<OneD, int>(pivotSize);
 
                             Lapack::Dgetrf(m, n, A.get(), m, m_ipivot.get(), info);
 
                             if( info < 0 )
                             {
                                 std::string message = "ERROR: The " + boost::lexical_cast<std::string>(-info) + "th parameter had an illegal parameter for dgetrf";
                                 ASSERTL0(false, message.c_str());
                             }
                             else if( info > 0 )
                             {
                                 std::string message = "ERROR: Element u_" + boost::lexical_cast<std::string>(info) +   boost::lexical_cast<std::string>(info) + " is 0 from dgetrf";
                                 ASSERTL0(false, message.c_str());
                             } 
                         }
                         break;
                     case eDIAGONAL:
                         for(unsigned int i = 0; i < theA.GetColumns(); ++i)
                         {
                             A[i] = 1.0/theA(i,i);
                         }
                         break;
                     case eUPPER_TRIANGULAR:
                     case eLOWER_TRIANGULAR:
                         break;
                     case eSYMMETRIC:
                         {
                             int info = 0;
                             int pivotSize = theA.GetRows();
                             m_ipivot = Array<OneD, int>(pivotSize);
                             
                             Lapack::Dsptrf('U', theA.GetRows(), A.get(), m_ipivot.get(), info);
                             
                             if( info < 0 )
                             {
                                 std::string message = "ERROR: The " + boost::lexical_cast<std::string>(-info) + "th parameter had an illegal parameter for dsptrf";
                                 ASSERTL0(false, message.c_str());
                             }
                             else if( info > 0 )
                             {
                                 std::string message = "ERROR: Element u_" + boost::lexical_cast<std::string>(info) +   boost::lexical_cast<std::string>(info) + " is 0 from dsptrf";
                                 ASSERTL0(false, message.c_str());
                             }
                         }
                         break;
                     case ePOSITIVE_DEFINITE_SYMMETRIC:
                         {
                             int info = 0;
                             Lapack::Dpptrf('U', theA.GetRows(), A.get(), info);
                           
                             if( info < 0 )
                             {
                                 std::string message = "ERROR: The " + boost::lexical_cast<std::string>(-info) + "th parameter had an illegal parameter for dpptrf";
                                 ASSERTL0(false, message.c_str());
                             }
                             else if( info > 0 )
                             {
                                 std::string message = "ERROR: The leading minor of order " + boost::lexical_cast<std::string>(info) +  " is not positive definite from dpptrf";
                                 ASSERTL0(false, message.c_str());
                             }
                         }
                         break;
                     case eBANDED:
                         {
                             int M = n;
                             int N = n;
                             int KL = m_numberOfSubDiagonals;
                             int KU = m_numberOfSuperDiagonals;
                             
                             // The array we pass in to dgbtrf must have enough space for KL
                             // subdiagonals and KL+KU superdiagonals (see lapack users guide,
                             // in the section discussing band storage.
                             unsigned int requiredStorageSize = BandedMatrixFuncs::
                                 GetRequiredStorageSize(n, n, KL, KL+KU);
                             
                             unsigned int rawRows = KL+KU+1;
                             A = Array<OneD, double>(requiredStorageSize);
 
                             // Put the extra elements up front.
                             for(unsigned int i = 0; i < theA.GetColumns(); ++i)
                             {
                                 std::copy(theA.GetRawPtr() + i*rawRows, theA.GetRawPtr() + (i+1)*rawRows,
                                     A.get() + (i+1)*KL + i*rawRows);
                             }
                                    
                             int info = 0;
                             int pivotSize = theA.GetRows();
                             m_ipivot = Array<OneD, int>(pivotSize);
                             
                             Lapack::Dgbtrf(M, N, KL, KU, A.get(), 2*KL+KU+1, m_ipivot.get(), info);
                           
                             if( info < 0 )
                             {
                                 std::string message = "ERROR: The " + boost::lexical_cast<std::string>(-info) + "th parameter had an illegal parameter for dgbtrf";
                                 ASSERTL0(false, message.c_str());
                             }
                             else if( info > 0 )
                             {
                                 std::string message = "ERROR: Element u_" + boost::lexical_cast<std::string>(info) +   boost::lexical_cast<std::string>(info) + " is 0 from dgbtrf";
                                 ASSERTL0(false, message.c_str());
                             }
                         }
                         break;
                     case ePOSITIVE_DEFINITE_SYMMETRIC_BANDED:
                         {
                             ASSERTL1(m_numberOfSuperDiagonals==m_numberOfSuperDiagonals,
                                      std::string("Number of sub- and superdiagonals should ") + 
                                      std::string("be equal for a symmetric banded matrix"));
 
                             int KU = m_numberOfSuperDiagonals;
                             int info = 0;
                             Lapack::Dpbtrf('U', theA.GetRows(), KU, A.get(), KU+1, info);
                           
                             if( info < 0 )
                             {
                                 std::string message = "ERROR: The " + boost::lexical_cast<std::string>(-info) + "th parameter had an illegal parameter for dpbtrf";
                                 ASSERTL0(false, message.c_str());
                             }
                             else if( info > 0 )
                             {
                                 std::string message = "ERROR: The leading minor of order " + boost::lexical_cast<std::string>(info) +  " is not positive definite from dpbtrf";
                                 ASSERTL0(false, message.c_str());
                             }
                         }
                         break;
                     case eSYMMETRIC_BANDED:
                         NEKERROR(ErrorUtil::efatal, "Unhandled matrix type");
                         break;
                     case eUPPER_TRIANGULAR_BANDED:
                         NEKERROR(ErrorUtil::efatal, "Unhandled matrix type");
                         break;
                     case eLOWER_TRIANGULAR_BANDED:
                         NEKERROR(ErrorUtil::efatal, "Unhandled matrix type");
                         break;
                         
                     default:
                         NEKERROR(ErrorUtil::efatal, "Unhandled matrix type");
                 }
             }

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

inline

Definition at line 438 of file NekLinSys.hpp.

References n.

438 { return n; }

Nektar::LinearSystem::n

unsigned int n

Definition: NekLinSys.hpp:604

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

Nektar::LinearSystem::n

unsigned int n

Definition: NekLinSys.hpp:604

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

inline

Definition at line 391 of file NekLinSys.hpp.

References swap().

             {
                 LinearSystem temp(rhs);
                 swap(temp);
                 return *this;
             }

template<typename VectorType >

RawType<VectorType>::type Nektar::LinearSystem::Solve ( const VectorType & b )

inline

Definition at line 403 of file NekLinSys.hpp.

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

Referenced by Nektar::LibUtilities::Invert(), and Nektar::StdRegions::Invert().

             {
                 typename RawType<VectorType>::type x(ConsistentObjectAccess<VectorType>::const_reference(b).GetRows());
                 LinearSystemSolver::Solve(ConsistentObjectAccess<VectorType>::const_reference(b), x, m_matrixType,
                     m_ipivot, n, A, m_transposeFlag, m_numberOfSubDiagonals, m_numberOfSuperDiagonals);
                 return x;
             }    

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 A, m_ipivot, m_matrixType, m_numberOfSubDiagonals, m_numberOfSuperDiagonals, m_transposeFlag, n, and Nektar::LinearSystemSolver::Solve().

             {
                 LinearSystemSolver::Solve(ConsistentObjectAccess<BType>::const_reference(b), 
                       ConsistentObjectAccess<XType>::reference(x), m_matrixType,
                       m_ipivot, n, A, m_transposeFlag, m_numberOfSubDiagonals, m_numberOfSuperDiagonals);
             }

template<typename VectorType >

RawType<VectorType>::type Nektar::LinearSystem::SolveTranspose ( const VectorType & b )

inline

Definition at line 421 of file NekLinSys.hpp.

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

Referenced by Nektar::LibUtilities::MakeQuadratureWeights(), and Nektar::LibUtilities::MakeTetWeights().

             {
                 typename RawType<VectorType>::type x(ConsistentObjectAccess<VectorType>::const_reference(b).GetRows());
                 LinearSystemSolver::SolveTranspose(ConsistentObjectAccess<VectorType>::const_reference(b), x, m_matrixType,
                     m_ipivot, n, A, m_transposeFlag, m_numberOfSubDiagonals, m_numberOfSuperDiagonals);
                 return x;
             }    

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 A, m_ipivot, m_matrixType, m_numberOfSubDiagonals, m_numberOfSuperDiagonals, m_transposeFlag, n, and Nektar::LinearSystemSolver::SolveTranspose().

             {
                 LinearSystemSolver::SolveTranspose(ConsistentObjectAccess<BType>::const_reference(b), 
                                ConsistentObjectAccess<XType>::reference(x), m_matrixType,
                                m_ipivot, n, A, m_transposeFlag, m_numberOfSubDiagonals, m_numberOfSuperDiagonals);
             }

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=().

             {
                 std::swap(n, rhs.n);
                 std::swap(A, rhs.A);
                 std::swap(m_ipivot,rhs.m_ipivot);
                 std::swap(m_numberOfSubDiagonals, rhs.m_numberOfSubDiagonals);
                 std::swap(m_numberOfSuperDiagonals, rhs.m_numberOfSuperDiagonals);
                 std::swap(m_matrixType, rhs.m_matrixType);
                 std::swap(m_transposeFlag, rhs.m_transposeFlag);
             }