Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator 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 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.

{}

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.

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

Definition at line 437 of file NekLinSys.hpp.

References n.

Referenced by Solve(), and SolveTranspose().

{ return n; }
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, and n.

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

{
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 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, and n.

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

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

Member Data Documentation

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

Definition at line 605 of file NekLinSys.hpp.

Referenced by FactorMatrix(), LinearSystem(), Solve(), SolveTranspose(), 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().

Generated on Sun Mar 15 2015 00:11:48 for Nektar++ by   doxygen 1.8.1.2