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)
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 376 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 380 of file NekLinSys.hpp.

        : 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().size());
            CopyArray(theA->GetPtr(), A);
        }
 
        FactorMatrix(*theA);
    }

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

        : 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().size());
            CopyArray(theA.GetPtr(), A);
        }
 
        FactorMatrix(theA);
    }

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

Member Function Documentation

FactorMatrix()

template<typename MatrixType >

void Nektar::LinearSystem::FactorMatrix ( const MatrixType &  theA )

inlineprivate

Definition at line 502 of file NekLinSys.hpp.

    {
        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 " + std::to_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_" + std::to_string(info) +
                        std::to_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 " + std::to_string(-info) +
                        "th parameter had an illegal parameter for dsptrf";
                    NEKERROR(ErrorUtil::efatal, message.c_str());
                }
                else if (info > 0)
                {
                    std::string message =
                        "ERROR: Element u_" + std::to_string(info) +
                        std::to_string(info) + " is 0 from dsptrf";
                    NEKERROR(ErrorUtil::efatal, 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 " + std::to_string(-info) +
                        "th parameter had an illegal parameter for dpptrf";
                    NEKERROR(ErrorUtil::efatal, message.c_str());
                }
                else if (info > 0)
                {
                    std::string message =
                        "ERROR: The leading minor of order " +
                        std::to_string(info) +
                        " is not positive definite from dpptrf";
                    NEKERROR(ErrorUtil::efatal, 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 " + std::to_string(-info) +
                        "th parameter had an illegal parameter for dgbtrf";
                    NEKERROR(ErrorUtil::efatal, message.c_str());
                }
                else if (info > 0)
                {
                    std::string message =
                        "ERROR: Element u_" + std::to_string(info) +
                        std::to_string(info) + " is 0 from dgbtrf";
                    NEKERROR(ErrorUtil::efatal, 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 " + std::to_string(-info) +
                        "th parameter had an illegal parameter for dpbtrf";
                    NEKERROR(ErrorUtil::efatal, message.c_str());
                }
                else if (info > 0)
                {
                    std::string message =
                        "ERROR: The leading minor of order " +
                        std::to_string(info) +
                        " is not positive definite from dpbtrf";
                    NEKERROR(ErrorUtil::efatal, 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");
        }
    }

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

    {
        return n;
    }

References n.

GetRows()

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

inline

Definition at line 492 of file NekLinSys.hpp.

    {
        return n;
    }

References n.

Referenced by Solve(), and SolveTranspose().

operator=()

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

inline

Definition at line 438 of file NekLinSys.hpp.

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

References swap().

Solve() [1/2]

template<typename BType , typename XType >

void Nektar::LinearSystem::Solve ( const BType &  b, XType &  x  ) const

inline

Definition at line 460 of file NekLinSys.hpp.

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

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

    {
        RawType_t<VectorType> 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;
    }

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

Referenced by Nektar::BandedMatrixVectorMultiplicationUnitTests::BOOST_AUTO_TEST_CASE(), Nektar::LinearSystemUnitTests::BOOST_AUTO_TEST_CASE(), and Nektar::TriangularMatrixVectorMultiplicationUnitTests::BOOST_AUTO_TEST_CASE().

SolveTranspose() [1/2]

template<typename BType , typename XType >

void Nektar::LinearSystem::SolveTranspose ( const BType &  b, XType &  x  ) const

inline

Definition at line 483 of file NekLinSys.hpp.

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

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

    {
        RawType_t<VectorType> 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;
    }

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

Referenced by Nektar::LinearSystemUnitTests::BOOST_AUTO_TEST_CASE(), Nektar::TriangularMatrixVectorMultiplicationUnitTests::BOOST_AUTO_TEST_CASE(), and Nektar::LibUtilities::NodalUtil::GetWeights().

swap()

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

inlineprivate

Definition at line 680 of file NekLinSys.hpp.

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

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

Referenced by swap().

m_ipivot

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

private

Definition at line 693 of file NekLinSys.hpp.

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

m_matrixType

MatrixStorage Nektar::LinearSystem::m_matrixType

private

Definition at line 696 of file NekLinSys.hpp.

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

m_numberOfSubDiagonals

unsigned int Nektar::LinearSystem::m_numberOfSubDiagonals

private

Definition at line 694 of file NekLinSys.hpp.

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

m_numberOfSuperDiagonals

unsigned int Nektar::LinearSystem::m_numberOfSuperDiagonals

private

Definition at line 695 of file NekLinSys.hpp.

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

m_transposeFlag

char Nektar::LinearSystem::m_transposeFlag

private

Definition at line 697 of file NekLinSys.hpp.

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

n

unsigned int Nektar::LinearSystem::n

private

Definition at line 691 of file NekLinSys.hpp.

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