A global linear system. More...

#include <GlobalLinSysIterative.h>

Inheritance diagram for Nektar::MultiRegions::GlobalLinSysIterative:

Public Member Functions
	GlobalLinSysIterative (const GlobalLinSysKey &pKey, const std::weak_ptr< ExpList > &pExpList, const std::shared_ptr< AssemblyMap > &pLocToGloMap)
	Constructor for full direct matrix solve.

	~GlobalLinSysIterative () override

void	DoMatrixMultiply (const Array< OneD, NekDouble > &pInput, Array< OneD, NekDouble > &pOutput)

void	DoPreconditionerFlag (const Array< OneD, NekDouble > &pInput, Array< OneD, NekDouble > &pOutput, const bool &isLocal=false)

Public Member Functions inherited from Nektar::MultiRegions::GlobalLinSys
	GlobalLinSys (const GlobalLinSysKey &pKey, const std::weak_ptr< ExpList > &pExpList, const std::shared_ptr< AssemblyMap > &pLocToGloMap)
	Constructor for full direct matrix solve.

virtual	~GlobalLinSys ()

const GlobalLinSysKey &	GetKey (void) const
	Returns the key associated with the system.

const std::weak_ptr< ExpList > &	GetLocMat (void) const

void	InitObject ()

void	Initialise (const std::shared_ptr< AssemblyMap > &pLocToGloMap)

void	Solve (const Array< OneD, const NekDouble > &in, Array< OneD, NekDouble > &out, const AssemblyMapSharedPtr &locToGloMap, const Array< OneD, const NekDouble > &dirForcing=NullNekDouble1DArray)
	Solve the linear system for given input and output vectors using a specified local to global map.

std::shared_ptr< GlobalLinSys >	GetSharedThisPtr ()
	Returns a shared pointer to the current object.

int	GetNumBlocks ()

DNekScalMatSharedPtr	GetBlock (unsigned int n)

void	DropBlock (unsigned int n)

DNekScalBlkMatSharedPtr	GetStaticCondBlock (unsigned int n)

void	DropStaticCondBlock (unsigned int n)

void	SolveLinearSystem (const int pNumRows, const Array< OneD, const NekDouble > &pInput, Array< OneD, NekDouble > &pOutput, const AssemblyMapSharedPtr &locToGloMap, const int pNumDir=0)
	Solve the linear system for given input and output vectors.

Protected Member Functions
void	DoProjection (const int pNumRows, const Array< OneD, const NekDouble > &pInput, Array< OneD, NekDouble > &pOutput, const int pNumDir, const bool isAconjugate)
	projection technique

bool	isNonSymmetricLinSys (StdRegions::MatrixType mt)

virtual void	v_UniqueMap ()=0

virtual void	v_DoMatrixMultiply (const Array< OneD, NekDouble > &pInput, Array< OneD, NekDouble > &pOutput)=0

Protected Member Functions inherited from Nektar::MultiRegions::GlobalLinSys
virtual void	v_Solve (const Array< OneD, const NekDouble > &in, Array< OneD, NekDouble > &out, const AssemblyMapSharedPtr &locToGloMap, const Array< OneD, const NekDouble > &dirForcing=NullNekDouble1DArray)=0
	Solve a linear system based on mapping.

virtual void	v_SolveLinearSystem (const int pNumRows, const Array< OneD, const NekDouble > &pInput, Array< OneD, NekDouble > &pOutput, const AssemblyMapSharedPtr &locToGloMap, const int pNumDir)=0
	Solve a basic matrix system.

virtual void	v_InitObject ()

virtual void	v_Initialise (const std::shared_ptr< AssemblyMap > &pLocToGloMap)

virtual int	v_GetNumBlocks ()
	Get the number of blocks in this system.

virtual DNekScalMatSharedPtr	v_GetBlock (unsigned int n)
	Retrieves the block matrix from n-th expansion using the matrix key provided by the m_linSysKey.

virtual void	v_DropBlock (unsigned int n)
	Releases the local block matrix from NekManager of n-th expansion using the matrix key provided by the m_linSysKey.

virtual DNekScalBlkMatSharedPtr	v_GetStaticCondBlock (unsigned int n)
	Retrieves a the static condensation block matrices from n-th expansion using the matrix key provided by the m_linSysKey.

virtual void	v_DropStaticCondBlock (unsigned int n)
	Releases the static condensation block matrices from NekManager of n-th expansion using the matrix key provided by the m_linSysKey.

PreconditionerSharedPtr	CreatePrecon (AssemblyMapSharedPtr asmMap)
	Create a preconditioner object from the parameters defined in the supplied assembly map.

Protected Attributes
Array< OneD, int >	m_map
	Global to universal unique map.

NekDouble	m_rhs_magnitude
	dot product of rhs to normalise stopping criterion

NekDouble	m_rhs_mag_sm
	cnt to how many times rhs_magnitude is called

PreconditionerSharedPtr	m_precon

std::string	m_precontype

int	m_totalIterations

bool	m_useProjection
	Whether to apply projection technique.

bool	m_isConjugateGradient

bool	m_isLocalSolver

bool	m_root
	Root if parallel.

std::string	m_linSysIterSolver
	Iterative solver: Conjugate Gradient, GMRES.

std::vector< Array< OneD, NekDouble > >	m_prevLinSol
	Storage for solutions to previous linear problems.

std::vector< Array< OneD, NekDouble > >	m_prevBasis

DNekMatSharedPtr	m_coeffMatrix

Array< OneD, NekDouble >	m_coeffMatrixFactor

Array< OneD, int >	m_ipivot

int	m_numSuccessiveRHS

bool	m_isAconjugate

std::string	m_matrixType

int	m_numPrevSols

bool	m_isAbsoluteTolerance

LibUtilities::NekSysOperators	m_NekSysOp

LibUtilities::NekLinSysIterSharedPtr	m_linsol

Protected Attributes inherited from Nektar::MultiRegions::GlobalLinSys
const GlobalLinSysKey	m_linSysKey
	Key associated with this linear system.

const std::weak_ptr< ExpList >	m_expList
	Local Matrix System.

const std::map< int, RobinBCInfoSharedPtr >	m_robinBCInfo
	Robin boundary info.

bool	m_verbose

Static Protected Attributes
static std::string	IteratSolverlookupIds []

static std::string	IteratSolverdef

Private Member Functions
void	UpdateKnownSolutions (const int pGlobalBndDofs, const Array< OneD, const NekDouble > &pSolution, const int pNumDirBndDofs, const bool isAconjugate)

int	ResetKnownSolutionsToLatestOne ()

void	DoMatrixMultiplyFlag (const Array< OneD, NekDouble > &pInput, Array< OneD, NekDouble > &pOutput, const bool &controlFlag)

void	DoAssembleLocFlag (const Array< OneD, NekDouble > &pInput, Array< OneD, NekDouble > &pOutput, const bool &ZeroDir)

Detailed Description

A global linear system.

Solves a linear system using iterative methods.

Definition at line 48 of file GlobalLinSysIterative.h.

Constructor & Destructor Documentation

◆ GlobalLinSysIterative()

Nektar::MultiRegions::GlobalLinSysIterative::GlobalLinSysIterative	(	const GlobalLinSysKey &	pKey,
		const std::weak_ptr< ExpList > &	pExpList,
		const std::shared_ptr< AssemblyMap > &	pLocToGloMap
	)

Constructor for full direct matrix solve.

Constructor for full iterative matrix solve.

Definition at line 67 of file GlobalLinSysIterative.cpp.

    : GlobalLinSys(pKey, pExpList, pLocToGloMap),
      m_rhs_magnitude(NekConstants::kNekUnsetDouble), m_rhs_mag_sm(0.9),
      m_precon(NullPreconditionerSharedPtr), m_totalIterations(0),
      m_useProjection(false), m_numPrevSols(0)
{
    m_isAbsoluteTolerance = pLocToGloMap->IsAbsoluteTolerance();
    m_linSysIterSolver    = pLocToGloMap->GetLinSysIterSolver();
 
    LibUtilities::CommSharedPtr vComm =
        m_expList.lock()->GetComm()->GetRowComm();
    m_root = (vComm->GetRank()) ? false : true;
 
    m_numSuccessiveRHS = pLocToGloMap->GetSuccessiveRHS();
    m_isAconjugate     = m_numSuccessiveRHS > 0;
    m_numSuccessiveRHS = std::abs(m_numSuccessiveRHS);
    m_useProjection    = m_numSuccessiveRHS > 0;
    m_isConjugateGradient =
        m_linSysIterSolver.find("ConjugateGradient") != std::string::npos;
    m_isLocalSolver = m_linSysIterSolver.find("Loc") != std::string::npos;
 
    // Check for nonsymmetric matrix and switch to GMRES, if not already used
    auto matrixType = pKey.GetMatrixType();
    if (isNonSymmetricLinSys(matrixType) && m_isConjugateGradient)
    {
        m_linSysIterSolver = m_isLocalSolver ? "GMRESLoc" : "GMRES";
 
        WARNINGL0(
            false,
            "Detected ConjugateGradient solver and a "
            "Advection-type matrix. "
            "Switchted to a GMRES solver for this non-symmetric matrix type. "
            "Change SolverInfo 'LinSysIterSolver' to GMRES in the session file "
            "to suppress this warning.");
    }
 
    if (m_isAconjugate && !m_linSysIterSolver.compare("GMRES") &&
        !m_linSysIterSolver.compare("GMRESLoc"))
    {
        WARNINGL0(false, "To use A-conjugate projection, the matrix "
                         "should be symmetric positive definite.");
    }
 
    m_NekSysOp.DefineNekSysLhsEval(&GlobalLinSysIterative::DoMatrixMultiplyFlag,
                                   this);
    m_NekSysOp.DefineNekSysPrecon(&GlobalLinSysIterative::DoPreconditionerFlag,
                                  this);
    m_NekSysOp.DefineAssembleLoc(&GlobalLinSysIterative::DoAssembleLocFlag,
                                 this);
}

References Nektar::LibUtilities::NekSysOperators::DefineAssembleLoc(), Nektar::LibUtilities::NekSysOperators::DefineNekSysLhsEval(), Nektar::LibUtilities::NekSysOperators::DefineNekSysPrecon(), DoAssembleLocFlag(), DoMatrixMultiplyFlag(), DoPreconditionerFlag(), Nektar::MultiRegions::GlobalMatrixKey::GetMatrixType(), isNonSymmetricLinSys(), Nektar::MultiRegions::GlobalLinSys::m_expList, m_isAbsoluteTolerance, m_isAconjugate, m_isConjugateGradient, m_isLocalSolver, m_linSysIterSolver, m_NekSysOp, m_numSuccessiveRHS, m_root, m_useProjection, and WARNINGL0.

◆ ~GlobalLinSysIterative()

Nektar::MultiRegions::GlobalLinSysIterative::~GlobalLinSysIterative ( )

override

Definition at line 120 of file GlobalLinSysIterative.cpp.

121{

122}

Member Function Documentation

◆ DoAssembleLocFlag()

void Nektar::MultiRegions::GlobalLinSysIterative::DoAssembleLocFlag	(	const Array< OneD, NekDouble > &	pInput,
		Array< OneD, NekDouble > &	pOutput,
		const bool &	ZeroDir
	)

inlineprivate

Definition at line 149 of file GlobalLinSysIterative.h.

    {
        m_precon->DoAssembleLoc(pInput, pOutput, ZeroDir);
    }

References m_precon.

Referenced by GlobalLinSysIterative().

◆ DoMatrixMultiply()

void Nektar::MultiRegions::GlobalLinSysIterative::DoMatrixMultiply	(	const Array< OneD, NekDouble > &	pInput,
		Array< OneD, NekDouble > &	pOutput
	)

inline

Definition at line 58 of file GlobalLinSysIterative.h.

    {
        v_DoMatrixMultiply(pInput, pOutput);
    }

References v_DoMatrixMultiply().

◆ DoMatrixMultiplyFlag()

void Nektar::MultiRegions::GlobalLinSysIterative::DoMatrixMultiplyFlag	(	const Array< OneD, NekDouble > &	pInput,
		Array< OneD, NekDouble > &	pOutput,
		const bool &	controlFlag
	)

inlineprivate

Definition at line 142 of file GlobalLinSysIterative.h.

    {
        v_DoMatrixMultiply(pInput, pOutput);
    }

References v_DoMatrixMultiply().

Referenced by DoProjection(), GlobalLinSysIterative(), and UpdateKnownSolutions().

◆ DoPreconditionerFlag()

void Nektar::MultiRegions::GlobalLinSysIterative::DoPreconditionerFlag	(	const Array< OneD, NekDouble > &	pInput,
		Array< OneD, NekDouble > &	pOutput,
		const bool &	isLocal = `false`
	)

inline

Definition at line 64 of file GlobalLinSysIterative.h.

    {
        m_precon->DoPreconditioner(pInput, pOutput, isLocal);
    }

References m_precon.

Referenced by GlobalLinSysIterative().

◆ DoProjection()

void Nektar::MultiRegions::GlobalLinSysIterative::DoProjection	(	const int	nGlobal,
		const Array< OneD, const NekDouble > &	pInput,
		Array< OneD, NekDouble > &	pOutput,
		const int	nDir,
		const bool	isAconjugate
	)

protected

projection technique

This method implements projection techniques in order to speed up successive linear solves with right-hand sides arising from time-dependent discretisations. (P.F.Fischer, Comput. Methods Appl. Mech. Engrg. 163, 1998)

Definition at line 130 of file GlobalLinSysIterative.cpp.

{
    int numIterations = 0;
    if (m_numPrevSols == 0)
    {
        // no previous solutions found
        numIterations = m_linsol->SolveSystem(nGlobal, pInput, pOutput, nDir);
    }
    else
    {
        // Get the communicator for performing data exchanges
        LibUtilities::CommSharedPtr vComm =
            m_expList.lock()->GetComm()->GetRowComm();
 
        // Get vector sizes
        int nNonDir = nGlobal - nDir;
 
        // check the input vector (rhs) is not zero
        Array<OneD, NekDouble> tmp;
 
        NekDouble rhsNorm =
            Vmath::Dot2(nNonDir, pInput + nDir, pInput + nDir, m_map + nDir);
 
        vComm->AllReduce(rhsNorm, Nektar::LibUtilities::ReduceSum);
 
        NekDouble tol = m_linsol->GetNekLinSysTolerance();
        if (rhsNorm < tol * tol * m_rhs_magnitude)
        {
            Vmath::Zero(nNonDir, tmp = pOutput + nDir, 1);
            if (m_verbose && m_root)
            {
                cout << "No iterations made"
                     << " using tolerance of " << tol
                     << " (error = " << sqrt(rhsNorm / m_rhs_magnitude)
                     << ", rhs_mag = " << sqrt(m_rhs_magnitude) << ")" << endl;
            }
            return;
        }
 
        // Create NekVector wrappers for linear algebra operations
        NekVector<NekDouble> b(nNonDir, pInput + nDir, eWrapper);
        NekVector<NekDouble> x(nNonDir, tmp = pOutput + nDir, eWrapper);
        // Allocate array storage
        Array<OneD, NekDouble> px_s(nGlobal, 0.0);
        Array<OneD, NekDouble> pb_s(nGlobal, 0.0);
        Array<OneD, NekDouble> tmpAx_s(nGlobal, 0.0);
        Array<OneD, NekDouble> tmpx_s(nGlobal, 0.0);
 
        NekVector<NekDouble> pb(nNonDir, tmp = pb_s + nDir, eWrapper);
        NekVector<NekDouble> px(nNonDir, tmp = px_s + nDir, eWrapper);
        NekVector<NekDouble> tmpAx(nNonDir, tmp = tmpAx_s + nDir, eWrapper);
        NekVector<NekDouble> tmpx(nNonDir, tmp = tmpx_s + nDir, eWrapper);
 
        // notation follows the paper cited:
        // \alpha_i = \tilda{x_i}^T b^n
        // projected x, px = \sum \alpha_i \tilda{x_i}
 
        Array<OneD, NekDouble> alpha(m_prevBasis.size(), 0.0);
        Array<OneD, NekDouble> alphaback(m_prevBasis.size(), 0.0);
        for (int i = 0; i < m_prevBasis.size(); i++)
        {
            alpha[i] = Vmath::Dot2(nNonDir, m_prevBasis[i], pInput + nDir,
                                   m_map + nDir);
        }
        vComm->AllReduce(alpha, Nektar::LibUtilities::ReduceSum);
        int n = m_prevBasis.size(), info = -1;
        Vmath::Vcopy(m_prevBasis.size(), alpha, 1, alphaback, 1);
        Lapack::Dsptrs('U', n, 1, m_coeffMatrixFactor.data(), m_ipivot.data(),
                       alpha.data(), n, info);
        if (info != 0)
        {
            // Dsptrs fails, only keep the latest solution
            int latest = ResetKnownSolutionsToLatestOne();
            alpha[0]   = alphaback[latest];
        }
        for (int i = 0; i < m_prevBasis.size(); ++i)
        {
            NekVector<NekDouble> xi(nNonDir, m_prevLinSol[i], eWrapper);
            px += alpha[i] * xi;
        }
 
        // pb = b^n - A px
        Vmath::Vcopy(nNonDir, pInput.data() + nDir, 1, pb_s.data() + nDir, 1);
 
        DoMatrixMultiplyFlag(px_s, tmpAx_s, false);
 
        pb -= tmpAx;
 
        if (m_verbose)
        {
            if (m_root)
            {
                cout << "SuccessiveRHS: " << m_prevBasis.size()
                     << "-bases projection reduces L2-norm of RHS from "
                     << std::sqrt(rhsNorm) << " to ";
            }
            NekDouble tmprhsNorm =
                Vmath::Dot2(nNonDir, pb_s + nDir, pb_s + nDir, m_map + nDir);
            vComm->AllReduce(tmprhsNorm, Nektar::LibUtilities::ReduceSum);
            if (m_root)
            {
                cout << std::sqrt(tmprhsNorm) << endl;
            }
        }
 
        // solve the system with projected rhs
        numIterations = m_linsol->SolveSystem(nGlobal, pb_s, tmpx_s, nDir);
 
        // remainder solution + projection of previous solutions
        x = tmpx + px;
    }
    // save the auxiliary solution to prev. known solutions
    if (numIterations)
    {
        UpdateKnownSolutions(nGlobal, pOutput, nDir, isAconjugate);
    }
}

References DoMatrixMultiplyFlag(), Vmath::Dot2(), Lapack::Dsptrs(), Nektar::eWrapper, m_coeffMatrixFactor, Nektar::MultiRegions::GlobalLinSys::m_expList, m_ipivot, m_linsol, m_map, m_numPrevSols, m_prevBasis, m_prevLinSol, m_rhs_magnitude, m_root, Nektar::MultiRegions::GlobalLinSys::m_verbose, Nektar::LibUtilities::ReduceSum, ResetKnownSolutionsToLatestOne(), tinysimd::sqrt(), UpdateKnownSolutions(), Vmath::Vcopy(), and Vmath::Zero().

Referenced by Nektar::MultiRegions::GlobalLinSysIterativeFull::v_SolveLinearSystem(), and Nektar::MultiRegions::GlobalLinSysIterativeStaticCond::v_SolveLinearSystem().

◆ isNonSymmetricLinSys()

bool Nektar::MultiRegions::GlobalLinSysIterative::isNonSymmetricLinSys ( StdRegions::MatrixType mt )

protected

Check if the matrix type is non-symmetric to allow swtiching to a capable iterative solver.

Parameters

mt	type of the Matrix using StdRegions::MatrixType

Returns: boolean if it is a non-symmetric type

Definition at line 488 of file GlobalLinSysIterative.cpp.

{
    constexpr std::array<StdRegions::MatrixType, 4> nonSymmetric = {
        StdRegions::eLinearAdvection,
        StdRegions::eLinearAdvectionReaction,
        StdRegions::eLinearAdvectionDiffusionReaction,
        StdRegions::eLinearAdvectionDiffusionReactionGJP,
    };
    return std::find(nonSymmetric.begin(), nonSymmetric.end(), mt) !=
           nonSymmetric.end();
}

References Nektar::StdRegions::eLinearAdvection, Nektar::StdRegions::eLinearAdvectionDiffusionReaction, Nektar::StdRegions::eLinearAdvectionDiffusionReactionGJP, and Nektar::StdRegions::eLinearAdvectionReaction.

Referenced by GlobalLinSysIterative().

◆ ResetKnownSolutionsToLatestOne()

int Nektar::MultiRegions::GlobalLinSysIterative::ResetKnownSolutionsToLatestOne ( )

private

Definition at line 250 of file GlobalLinSysIterative.cpp.

{
    if (m_numPrevSols == 0)
    {
        return -1;
    }
    int latest = (m_numPrevSols - 1 + m_numSuccessiveRHS) % m_numSuccessiveRHS;
    Array<OneD, NekDouble> b = m_prevBasis[latest];
    Array<OneD, NekDouble> x = m_prevLinSol[latest];
    m_prevBasis.clear();
    m_prevLinSol.clear();
    m_prevBasis.push_back(b);
    m_prevLinSol.push_back(x);
    m_numPrevSols = 1;
    return latest;
}

References m_numPrevSols, m_numSuccessiveRHS, m_prevBasis, and m_prevLinSol.

Referenced by DoProjection(), and UpdateKnownSolutions().

◆ UpdateKnownSolutions()

void Nektar::MultiRegions::GlobalLinSysIterative::UpdateKnownSolutions	(	const int	nGlobal,
		const Array< OneD, const NekDouble > &	newX,
		const int	nDir,
		const bool	isAconjugate
	)

private

Updates the storage of previously known solutions. Performs normalisation of input vector wrt A-norm.

Definition at line 271 of file GlobalLinSysIterative.cpp.

{
    // Get vector sizes
    int nNonDir        = nGlobal - nDir;
    int insertLocation = m_numPrevSols % m_numSuccessiveRHS;
    int fullbuffer     = (m_prevBasis.size() == m_numSuccessiveRHS);
 
    // Get the communicator for performing data exchanges
    LibUtilities::CommSharedPtr vComm =
        m_expList.lock()->GetComm()->GetRowComm();
 
    Array<OneD, NekDouble> tmpAx_s(nGlobal, 0.0);
    Array<OneD, NekDouble> y_s(m_prevBasis.size() - fullbuffer, 0.0);
    Array<OneD, NekDouble> invMy_s(y_s.size(), 0.0);
    Array<OneD, int> ipivot(m_numSuccessiveRHS);
    Array<OneD, NekDouble> tmp, newBasis;
 
    DoMatrixMultiplyFlag(newX, tmpAx_s, false);
 
    if (isAconjugate)
    {
        newBasis = newX + nDir;
    }
    else
    {
        newBasis = tmpAx_s + nDir;
    }
 
    // Check the solution is non-zero
    NekDouble solNorm =
        Vmath::Dot2(nNonDir, newBasis, tmpAx_s + nDir, m_map + nDir);
    vComm->AllReduce(solNorm, Nektar::LibUtilities::ReduceSum);
 
    if (solNorm < NekConstants::kNekSparseNonZeroTol)
    {
        return;
    }
 
    // normalisation of A x
    Vmath::Smul(nNonDir, 1.0 / sqrt(solNorm), tmpAx_s + nDir, 1,
                tmp = tmpAx_s + nDir, 1);
 
    for (int i = 0; i < m_prevBasis.size(); ++i)
    {
        if (i == insertLocation)
        {
            continue;
        }
        int skip = i > insertLocation;
        y_s[i - skip] =
            Vmath::Dot2(nNonDir, m_prevBasis[i], tmpAx_s + nDir, m_map + nDir);
    }
    vComm->AllReduce(y_s, Nektar::LibUtilities::ReduceSum);
 
    // check if linearly dependent
    DNekMatSharedPtr tilCoeffMatrix;
    if (fullbuffer && m_numSuccessiveRHS > 1)
    {
        tilCoeffMatrix = MemoryManager<DNekMat>::AllocateSharedPtr(
            m_numSuccessiveRHS - 1, m_numSuccessiveRHS - 1, 0.0, eSYMMETRIC);
        for (int i = 0; i < m_numSuccessiveRHS; ++i)
        {
            if (i == insertLocation)
            {
                continue;
            }
            int iskip = i > insertLocation;
            for (int j = i; j < m_numSuccessiveRHS; ++j)
            {
                if (j == insertLocation)
                {
                    continue;
                }
                int jskip = j > insertLocation;
                tilCoeffMatrix->SetValue(i - iskip, j - jskip,
                                         m_coeffMatrix->GetValue(i, j));
            }
        }
    }
    else if (!fullbuffer && m_prevBasis.size())
    {
        tilCoeffMatrix = MemoryManager<DNekMat>::AllocateSharedPtr(
            m_prevBasis.size(), m_prevBasis.size(), 0.0, eSYMMETRIC);
        Vmath::Vcopy(tilCoeffMatrix->GetStorageSize(), m_coeffMatrix->GetPtr(),
                     1, tilCoeffMatrix->GetPtr(), 1);
    }
 
    int n, info1 = 0, info2 = 0, info3 = 0;
    if (y_s.size())
    {
        n = tilCoeffMatrix->GetRows();
        Array<OneD, NekDouble> tilCoeffMatrixFactor(
            tilCoeffMatrix->GetStorageSize());
        Vmath::Vcopy(tilCoeffMatrix->GetStorageSize(), tilCoeffMatrix->GetPtr(),
                     1, tilCoeffMatrixFactor, 1);
        Lapack::Dsptrf('U', n, tilCoeffMatrixFactor.data(), ipivot.data(),
                       info1);
        if (info1 == 0)
        {
            Vmath::Vcopy(n, y_s, 1, invMy_s, 1);
            Lapack::Dsptrs('U', n, 1, tilCoeffMatrixFactor.data(),
                           ipivot.data(), invMy_s.data(), n, info2);
        }
    }
    if (info1 || info2)
    {
        int latest     = ResetKnownSolutionsToLatestOne();
        y_s[0]         = y_s[latest - (latest > insertLocation)];
        invMy_s[0]     = y_s[0];
        insertLocation = m_numPrevSols % m_numSuccessiveRHS;
        fullbuffer     = (m_prevBasis.size() == m_numSuccessiveRHS);
    }
    NekDouble residual = 1.;
    NekDouble epsilon  = 10. * NekConstants::kNekZeroTol;
    for (int i = 0; i < m_prevBasis.size() - fullbuffer; i++)
    {
        residual -= y_s[i] * invMy_s[i];
    }
    if (m_verbose && m_root)
    {
        cout << "SuccessiveRHS: residual " << residual;
    }
    if (residual < epsilon)
    {
        if (m_verbose && m_root)
        {
            cout << " < " << epsilon << ", reject" << endl;
        }
        return;
    }
 
    // calculate new coefficient matrix and its factor
    DNekMatSharedPtr newCoeffMatrix;
    if (fullbuffer)
    {
        newCoeffMatrix = MemoryManager<DNekMat>::AllocateSharedPtr(
            m_numSuccessiveRHS, m_numSuccessiveRHS, 0.0, eSYMMETRIC);
        Vmath::Vcopy(m_coeffMatrix->GetStorageSize(), m_coeffMatrix->GetPtr(),
                     1, newCoeffMatrix->GetPtr(), 1);
        newCoeffMatrix->SetValue(insertLocation, insertLocation, 1.);
        for (int i = 0; i < m_numSuccessiveRHS; ++i)
        {
            if (i == insertLocation)
            {
                continue;
            }
            int iskip = i > insertLocation;
            newCoeffMatrix->SetValue(insertLocation, i, y_s[i - iskip]);
        }
    }
    else
    {
        newCoeffMatrix = MemoryManager<DNekMat>::AllocateSharedPtr(
            m_prevBasis.size() + 1, m_prevBasis.size() + 1, 0.0, eSYMMETRIC);
        newCoeffMatrix->SetValue(insertLocation, insertLocation, 1.);
        for (int i = 0; i < m_prevBasis.size(); ++i)
        {
            newCoeffMatrix->SetValue(insertLocation, i, y_s[i]);
            for (int j = i; j < m_prevBasis.size(); ++j)
            {
                newCoeffMatrix->SetValue(i, j, m_coeffMatrix->GetValue(i, j));
            }
        }
    }
    n = newCoeffMatrix->GetRows();
    Array<OneD, NekDouble> coeffMatrixFactor(newCoeffMatrix->GetStorageSize());
    Vmath::Vcopy(newCoeffMatrix->GetStorageSize(), newCoeffMatrix->GetPtr(), 1,
                 coeffMatrixFactor, 1);
    Lapack::Dsptrf('U', n, coeffMatrixFactor.data(), ipivot.data(), info3);
    if (info3)
    {
        if (m_verbose && m_root)
        {
            cout << " >= " << epsilon << ", reject (Dsptrf fails)" << endl;
        }
        return;
    }
    if (m_verbose && m_root)
    {
        cout << " >= " << epsilon << ", accept" << endl;
    }
 
    // if success, update basis, rhs, coefficient matrix, and its factor
    if (m_prevBasis.size() < m_numSuccessiveRHS)
    {
        m_prevBasis.push_back(tmp = tmpAx_s + nDir);
        if (isAconjugate)
        {
            m_prevLinSol.push_back(tmp);
        }
        else
        {
            Array<OneD, NekDouble> solution(nNonDir, 0.0);
            m_prevLinSol.push_back(solution);
        }
    }
    Vmath::Smul(nNonDir, 1. / sqrt(solNorm), tmp = newX + nDir, 1,
                m_prevLinSol[insertLocation], 1);
    if (!isAconjugate)
    {
        m_prevBasis[insertLocation] = tmpAx_s + nDir;
    }
    m_coeffMatrix       = newCoeffMatrix;
    m_coeffMatrixFactor = coeffMatrixFactor;
    m_ipivot            = ipivot;
    ++m_numPrevSols;
}

References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), DoMatrixMultiplyFlag(), Vmath::Dot2(), Lapack::Dsptrf(), Lapack::Dsptrs(), Nektar::eSYMMETRIC, Nektar::NekConstants::kNekSparseNonZeroTol, Nektar::NekConstants::kNekZeroTol, m_coeffMatrix, m_coeffMatrixFactor, Nektar::MultiRegions::GlobalLinSys::m_expList, m_ipivot, m_map, m_numPrevSols, m_numSuccessiveRHS, m_prevBasis, m_prevLinSol, m_root, Nektar::MultiRegions::GlobalLinSys::m_verbose, Nektar::LibUtilities::ReduceSum, ResetKnownSolutionsToLatestOne(), Vmath::Smul(), tinysimd::sqrt(), and Vmath::Vcopy().

Referenced by DoProjection().

◆ v_DoMatrixMultiply()

virtual void Nektar::MultiRegions::GlobalLinSysIterative::v_DoMatrixMultiply	(	const Array< OneD, NekDouble > &	pInput,
		Array< OneD, NekDouble > &	pOutput
	)

protectedpure virtual

Implemented in Nektar::MultiRegions::GlobalLinSysIterativeFull, and Nektar::MultiRegions::GlobalLinSysIterativeStaticCond.

Referenced by DoMatrixMultiply(), and DoMatrixMultiplyFlag().

◆ v_UniqueMap()

virtual void Nektar::MultiRegions::GlobalLinSysIterative::v_UniqueMap ( )

protectedpure virtual

Implemented in Nektar::MultiRegions::GlobalLinSysIterativeFull, and Nektar::MultiRegions::GlobalLinSysIterativeStaticCond.

Member Data Documentation

◆ IteratSolverdef

std::string Nektar::MultiRegions::GlobalLinSysIterative::IteratSolverdef

staticprotected

Initial value:

=
    LibUtilities::SessionReader::RegisterDefaultSolverInfo("LinSysIterSolver",
                                                           "ConjugateGradient")

Definition at line 119 of file GlobalLinSysIterative.h.

◆ IteratSolverlookupIds

std::string Nektar::MultiRegions::GlobalLinSysIterative::IteratSolverlookupIds

staticprotected

Initial value:

= {
    LibUtilities::SessionReader::RegisterEnumValue(
        "LinSysIterSolver", "ConjugateGradient",
        MultiRegions::eConjugateGradient),
    LibUtilities::SessionReader::RegisterEnumValue(
        "LinSysIterSolver", "ConjugateGradientLoc",
        MultiRegions::eConjugateGradient),
    LibUtilities::SessionReader::RegisterEnumValue("LinSysIterSolver", "GMRES",
                                                   MultiRegions::eGMRES),
    LibUtilities::SessionReader::RegisterEnumValue(
        "LinSysIterSolver", "GMRESLoc", MultiRegions::eGMRESLoc),
    LibUtilities::SessionReader::RegisterEnumValue(
        "LinSysIterSolver", "EvsDirect", MultiRegions::eDirectEigenValues),
}

Definition at line 118 of file GlobalLinSysIterative.h.

◆ m_coeffMatrix

DNekMatSharedPtr Nektar::MultiRegions::GlobalLinSysIterative::m_coeffMatrix

protected

Definition at line 105 of file GlobalLinSysIterative.h.

Referenced by UpdateKnownSolutions().

◆ m_coeffMatrixFactor

Array<OneD, NekDouble> Nektar::MultiRegions::GlobalLinSysIterative::m_coeffMatrixFactor

protected

Definition at line 106 of file GlobalLinSysIterative.h.

Referenced by DoProjection(), and UpdateKnownSolutions().

◆ m_ipivot

Array<OneD, int> Nektar::MultiRegions::GlobalLinSysIterative::m_ipivot

protected

Definition at line 107 of file GlobalLinSysIterative.h.

Referenced by DoProjection(), and UpdateKnownSolutions().

◆ m_isAbsoluteTolerance

bool Nektar::MultiRegions::GlobalLinSysIterative::m_isAbsoluteTolerance

protected

Definition at line 112 of file GlobalLinSysIterative.h.

Referenced by GlobalLinSysIterative(), Nektar::MultiRegions::GlobalLinSysIterativeStaticCond::v_PreSolve(), and Nektar::MultiRegions::GlobalLinSysIterativeFull::v_SolveLinearSystem().

◆ m_isAconjugate

bool Nektar::MultiRegions::GlobalLinSysIterative::m_isAconjugate

protected

Definition at line 109 of file GlobalLinSysIterative.h.

Referenced by GlobalLinSysIterative(), Nektar::MultiRegions::GlobalLinSysIterativeFull::v_SolveLinearSystem(), and Nektar::MultiRegions::GlobalLinSysIterativeStaticCond::v_SolveLinearSystem().

◆ m_isConjugateGradient

bool Nektar::MultiRegions::GlobalLinSysIterative::m_isConjugateGradient

protected

Definition at line 91 of file GlobalLinSysIterative.h.

Referenced by GlobalLinSysIterative().

◆ m_isLocalSolver

bool Nektar::MultiRegions::GlobalLinSysIterative::m_isLocalSolver

protected

Definition at line 94 of file GlobalLinSysIterative.h.

Referenced by GlobalLinSysIterative().

◆ m_linsol

LibUtilities::NekLinSysIterSharedPtr Nektar::MultiRegions::GlobalLinSysIterative::m_linsol

protected

Definition at line 116 of file GlobalLinSysIterative.h.

Referenced by DoProjection(), Nektar::MultiRegions::GlobalLinSysIterativeFull::Initialise(), Nektar::MultiRegions::GlobalLinSysIterativeFull::v_DoMatrixMultiply(), Nektar::MultiRegions::GlobalLinSysIterativeStaticCond::v_DoMatrixMultiply(), Nektar::MultiRegions::GlobalLinSysIterativeFull::v_SolveLinearSystem(), and Nektar::MultiRegions::GlobalLinSysIterativeStaticCond::v_SolveLinearSystem().

◆ m_linSysIterSolver

std::string Nektar::MultiRegions::GlobalLinSysIterative::m_linSysIterSolver

protected

Iterative solver: Conjugate Gradient, GMRES.

Definition at line 100 of file GlobalLinSysIterative.h.

Referenced by GlobalLinSysIterative(), Nektar::MultiRegions::GlobalLinSysIterativeFull::Initialise(), and Nektar::MultiRegions::GlobalLinSysIterativeStaticCond::v_SolveLinearSystem().

◆ m_map

Array<OneD, int> Nektar::MultiRegions::GlobalLinSysIterative::m_map

protected

Global to universal unique map.

Definition at line 73 of file GlobalLinSysIterative.h.

Referenced by DoProjection(), Nektar::MultiRegions::GlobalLinSysIterativeFull::Initialise(), UpdateKnownSolutions(), Nektar::MultiRegions::GlobalLinSysIterativeStaticCond::v_PreSolve(), Nektar::MultiRegions::GlobalLinSysIterativeStaticCond::v_SolveLinearSystem(), Nektar::MultiRegions::GlobalLinSysIterativeFull::v_UniqueMap(), and Nektar::MultiRegions::GlobalLinSysIterativeStaticCond::v_UniqueMap().

◆ m_matrixType

std::string Nektar::MultiRegions::GlobalLinSysIterative::m_matrixType

protected

Definition at line 110 of file GlobalLinSysIterative.h.

◆ m_NekSysOp

LibUtilities::NekSysOperators Nektar::MultiRegions::GlobalLinSysIterative::m_NekSysOp

protected

Definition at line 114 of file GlobalLinSysIterative.h.

Referenced by GlobalLinSysIterative(), Nektar::MultiRegions::GlobalLinSysIterativeFull::Initialise(), and Nektar::MultiRegions::GlobalLinSysIterativeStaticCond::v_SolveLinearSystem().

◆ m_numPrevSols

int Nektar::MultiRegions::GlobalLinSysIterative::m_numPrevSols

protected

Definition at line 111 of file GlobalLinSysIterative.h.

Referenced by DoProjection(), ResetKnownSolutionsToLatestOne(), and UpdateKnownSolutions().

◆ m_numSuccessiveRHS

int Nektar::MultiRegions::GlobalLinSysIterative::m_numSuccessiveRHS

protected

Definition at line 108 of file GlobalLinSysIterative.h.

Referenced by GlobalLinSysIterative(), ResetKnownSolutionsToLatestOne(), and UpdateKnownSolutions().

◆ m_precon

PreconditionerSharedPtr Nektar::MultiRegions::GlobalLinSysIterative::m_precon

protected

Definition at line 81 of file GlobalLinSysIterative.h.

◆ m_precontype

std::string Nektar::MultiRegions::GlobalLinSysIterative::m_precontype

protected

Definition at line 83 of file GlobalLinSysIterative.h.

◆ m_prevBasis

std::vector<Array<OneD, NekDouble> > Nektar::MultiRegions::GlobalLinSysIterative::m_prevBasis

protected

Definition at line 104 of file GlobalLinSysIterative.h.

Referenced by DoProjection(), ResetKnownSolutionsToLatestOne(), and UpdateKnownSolutions().

◆ m_prevLinSol

std::vector<Array<OneD, NekDouble> > Nektar::MultiRegions::GlobalLinSysIterative::m_prevLinSol

protected

Storage for solutions to previous linear problems.

Definition at line 103 of file GlobalLinSysIterative.h.

Referenced by DoProjection(), ResetKnownSolutionsToLatestOne(), and UpdateKnownSolutions().

◆ m_rhs_mag_sm

NekDouble Nektar::MultiRegions::GlobalLinSysIterative::m_rhs_mag_sm

protected

cnt to how many times rhs_magnitude is called

Definition at line 79 of file GlobalLinSysIterative.h.

Referenced by Nektar::MultiRegions::GlobalLinSysIterativeStaticCond::v_PreSolve().

◆ m_rhs_magnitude

NekDouble Nektar::MultiRegions::GlobalLinSysIterative::m_rhs_magnitude

protected

dot product of rhs to normalise stopping criterion

Definition at line 76 of file GlobalLinSysIterative.h.

Referenced by DoProjection(), Nektar::MultiRegions::GlobalLinSysIterativeStaticCond::v_PreSolve(), Nektar::MultiRegions::GlobalLinSysIterativeFull::v_SolveLinearSystem(), and Nektar::MultiRegions::GlobalLinSysIterativeStaticCond::v_SolveLinearSystem().

◆ m_root

bool Nektar::MultiRegions::GlobalLinSysIterative::m_root

protected

Root if parallel.

Definition at line 97 of file GlobalLinSysIterative.h.

Referenced by DoProjection(), GlobalLinSysIterative(), and UpdateKnownSolutions().

◆ m_totalIterations

int Nektar::MultiRegions::GlobalLinSysIterative::m_totalIterations

protected

Definition at line 85 of file GlobalLinSysIterative.h.

◆ m_useProjection

bool Nektar::MultiRegions::GlobalLinSysIterative::m_useProjection

protected

Whether to apply projection technique.

Definition at line 88 of file GlobalLinSysIterative.h.

Referenced by GlobalLinSysIterative(), Nektar::MultiRegions::GlobalLinSysIterativeFull::v_SolveLinearSystem(), and Nektar::MultiRegions::GlobalLinSysIterativeStaticCond::v_SolveLinearSystem().

Public Member Functions

Protected Member Functions

Protected Attributes

Static Protected Attributes

Private Member Functions

Detailed Description

Constructor & Destructor Documentation

◆ GlobalLinSysIterative()

◆ ~GlobalLinSysIterative()

Member Function Documentation

◆ DoAssembleLocFlag()

◆ DoMatrixMultiply()

◆ DoMatrixMultiplyFlag()

◆ DoPreconditionerFlag()

◆ DoProjection()

◆ isNonSymmetricLinSys()

◆ ResetKnownSolutionsToLatestOne()

◆ UpdateKnownSolutions()

◆ v_DoMatrixMultiply()

◆ v_UniqueMap()

Member Data Documentation

◆ IteratSolverdef

◆ IteratSolverlookupIds

◆ m_coeffMatrix

◆ m_coeffMatrixFactor

◆ m_ipivot

◆ m_isAbsoluteTolerance

◆ m_isAconjugate

◆ m_isConjugateGradient

◆ m_isLocalSolver

◆ m_linsol

◆ m_linSysIterSolver

◆ m_map

◆ m_matrixType

◆ m_NekSysOp

◆ m_numPrevSols

◆ m_numSuccessiveRHS

◆ m_precon

◆ m_precontype

◆ m_prevBasis

◆ m_prevLinSol

◆ m_rhs_mag_sm

◆ m_rhs_magnitude

◆ m_root

◆ m_totalIterations

◆ m_useProjection