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. More...

virtual	~GlobalLinSysIterative ()

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. More...

virtual	~GlobalLinSys ()

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

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. More...

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

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. More...

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

void	Set_Rhs_Magnitude (const NekVector< NekDouble > &pIn)

virtual void	v_UniqueMap ()=0

virtual void	v_SolveLinearSystem (const int pNumRows, const Array< OneD, const NekDouble > &pInput, Array< OneD, NekDouble > &pOutput, const AssemblyMapSharedPtr &locToGloMap, const int pNumDir) override
	Solve the matrix system. More...

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. More...

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. More...

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

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. More...

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. More...

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. More...

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

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

int	m_maxiter
	maximum iterations More...

NekDouble	m_tolerance
	Tolerance of iterative solver. More...

NekDouble	m_rhs_magnitude
	dot product of rhs to normalise stopping criterion More...

NekDouble	m_rhs_mag_sm
	cnt to how many times rhs_magnitude is called More...

PreconditionerSharedPtr	m_precon

MultiRegions::PreconditionerType	m_precontype

int	m_totalIterations

bool	m_useProjection
	Whether to apply projection technique. More...

bool	m_root
	Root if parallel. More...

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

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

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

int	m_numPrevSols

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. More...

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

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

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	DoPreconditionerFlag (const Array< OneD, NekDouble > &pInput, Array< OneD, NekDouble > &pOutput, const bool &controlFlag)

Detailed Description

A global linear system.

Solves a linear system using iterative methods.

Definition at line 50 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.

Definition at line 62 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_tolerance        = pLocToGloMap->GetIterativeTolerance();
     m_maxiter          = pLocToGloMap->GetMaxIterations();
     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;
  
     if (m_isAconjugate && 0 == m_linSysIterSolver.compare("GMRES"))
     {
         WARNINGL0(false, "To use A-conjugate projection, the matrix should be "
                          "symmetric positive definite.");
     }
 }

References tinysimd::abs(), Nektar::MultiRegions::GlobalLinSys::m_expList, m_isAconjugate, m_linSysIterSolver, m_maxiter, m_numSuccessiveRHS, m_root, m_tolerance, m_useProjection, and WARNINGL0.

◆ ~GlobalLinSysIterative()

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

virtual

Definition at line 90 of file GlobalLinSysIterative.cpp.

91 {

92 }

Member Function Documentation

◆ DoMatrixMultiplyFlag()

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

inlineprivate

Definition at line 136 of file GlobalLinSysIterative.h.

     {
         boost::ignore_unused(controlFlag);
  
         v_DoMatrixMultiply(pInput, pOutput);
     }

References v_DoMatrixMultiply().

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

◆ DoPreconditionerFlag()

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

inlineprivate

Definition at line 145 of file GlobalLinSysIterative.h.

     {
         boost::ignore_unused(controlFlag);
  
         m_precon->DoPreconditioner(pInput, pOutput);
     }

References m_precon.

Referenced by v_SolveLinearSystem().

◆ DoProjection()

void Nektar::MultiRegions::GlobalLinSysIterative::DoProjection	(	const int	nGlobal,
		const Array< OneD, const NekDouble > &	pInput,
		Array< OneD, NekDouble > &	pOutput,
		const int	nDir,
		const NekDouble	tol,
		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 151 of file GlobalLinSysIterative.cpp.

 {
     int numIterations = 0;
     if (0 == m_numPrevSols)
     {
         // no previous solutions found
         numIterations =
             m_linsol->SolveSystem(nGlobal, pInput, pOutput, nDir, tol);
     }
     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);
  
         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.get(), m_ipivot.get(),
                        alpha.get(), 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.get() + nDir, 1, pb_s.get() + 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, tol);
  
         // 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 v_SolveLinearSystem().

◆ ResetKnownSolutionsToLatestOne()

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

private

Definition at line 268 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().

◆ Set_Rhs_Magnitude()

void Nektar::MultiRegions::GlobalLinSysIterative::Set_Rhs_Magnitude ( const NekVector< NekDouble > & pIn )

protected

Definition at line 482 of file GlobalLinSysIterative.cpp.

 {
     Array<OneD, NekDouble> vExchange(1, 0.0);
     if (m_map.size() > 0)
     {
         vExchange[0] =
             Vmath::Dot2(pIn.GetDimension(), &pIn[0], &pIn[0], &m_map[0]);
     }
  
     m_expList.lock()->GetComm()->GetRowComm()->AllReduce(
         vExchange, Nektar::LibUtilities::ReduceSum);
  
     // To ensure that very different rhs values are not being
     // used in subsequent solvers such as the velocit solve in
     // INC NS. If this works we then need to work out a better
     // way to control this.
     NekDouble new_rhs_mag = (vExchange[0] > 1e-6) ? vExchange[0] : 1.0;
  
     if (m_rhs_magnitude == NekConstants::kNekUnsetDouble)
     {
         m_rhs_magnitude = new_rhs_mag;
     }
     else
     {
         m_rhs_magnitude = (m_rhs_mag_sm * (m_rhs_magnitude) +
                            (1.0 - m_rhs_mag_sm) * new_rhs_mag);
     }
 }

References Vmath::Dot2(), Nektar::NekVector< DataType >::GetDimension(), Nektar::NekConstants::kNekUnsetDouble, Nektar::MultiRegions::GlobalLinSys::m_expList, m_map, m_rhs_mag_sm, m_rhs_magnitude, and Nektar::LibUtilities::ReduceSum.

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

◆ 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 289 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 < 22.2 * 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.get(), ipivot.get(), info1);
         if (info1 == 0)
         {
             Vmath::Vcopy(n, y_s, 1, invMy_s, 1);
             Lapack::Dsptrs('U', n, 1, tilCoeffMatrixFactor.get(), ipivot.get(),
                            invMy_s.get(), 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.get(), ipivot.get(), 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::GlobalLinSysIterativeStaticCond, and Nektar::MultiRegions::GlobalLinSysIterativeFull.

Referenced by DoMatrixMultiplyFlag().

◆ v_SolveLinearSystem()

void Nektar::MultiRegions::GlobalLinSysIterative::v_SolveLinearSystem	(	const int	pNumRows,
		const Array< OneD, const NekDouble > &	pInput,
		Array< OneD, NekDouble > &	pOutput,
		const AssemblyMapSharedPtr &	locToGloMap,
		const int	pNumDir
	)

overrideprotectedvirtual

Solve the matrix system.

Implements Nektar::MultiRegions::GlobalLinSys.

Definition at line 97 of file GlobalLinSysIterative.cpp.

 {
     if (!m_linsol)
     {
         LibUtilities::CommSharedPtr v_Comm =
             m_expList.lock()->GetComm()->GetRowComm();
         LibUtilities::SessionReaderSharedPtr pSession =
             m_expList.lock()->GetSession();
  
         // Check such a module exists for this equation.
         ASSERTL0(LibUtilities::GetNekLinSysIterFactory().ModuleExists(
                      m_linSysIterSolver),
                  "NekLinSysIter '" + m_linSysIterSolver +
                      "' is not defined.\n");
         m_linsol = LibUtilities::GetNekLinSysIterFactory().CreateInstance(
             m_linSysIterSolver, pSession, v_Comm, nGlobal - nDir,
             LibUtilities::NekSysKey());
  
         m_NekSysOp.DefineNekSysLhsEval(
             &GlobalLinSysIterative::DoMatrixMultiplyFlag, this);
         m_NekSysOp.DefineNekSysPrecon(
             &GlobalLinSysIterative::DoPreconditionerFlag, this);
         m_linsol->SetSysOperators(m_NekSysOp);
         v_UniqueMap();
         m_linsol->setUniversalUniqueMap(m_map);
     }
  
     if (!m_precon)
     {
         m_precon = CreatePrecon(plocToGloMap);
         m_precon->BuildPreconditioner();
     }
  
     m_linsol->setRhsMagnitude(m_rhs_magnitude);
     if (m_useProjection)
     {
         DoProjection(nGlobal, pInput, pOutput, nDir, m_tolerance,
                      m_isAconjugate);
     }
     else
     {
         m_linsol->SolveSystem(nGlobal, pInput, pOutput, nDir, m_tolerance);
     }
 }

References ASSERTL0, Nektar::LibUtilities::NekFactory< tKey, tBase, tParam >::CreateInstance(), Nektar::MultiRegions::GlobalLinSys::CreatePrecon(), Nektar::LibUtilities::NekSysOperators::DefineNekSysLhsEval(), Nektar::LibUtilities::NekSysOperators::DefineNekSysPrecon(), DoMatrixMultiplyFlag(), DoPreconditionerFlag(), DoProjection(), Nektar::LibUtilities::GetNekLinSysIterFactory(), Nektar::MultiRegions::GlobalLinSys::m_expList, m_isAconjugate, m_linsol, m_linSysIterSolver, m_map, m_NekSysOp, m_precon, m_rhs_magnitude, m_tolerance, m_useProjection, and v_UniqueMap().

◆ v_UniqueMap()

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

protectedpure virtual

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

Referenced by v_SolveLinearSystem().

Member Data Documentation

◆ IteratSolverdef

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

staticprotected

Initial value:

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

Definition at line 106 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", "GMRES",
                                                   MultiRegions::eGMRES),
}

Definition at line 105 of file GlobalLinSysIterative.h.

◆ m_coeffMatrix

DNekMatSharedPtr Nektar::MultiRegions::GlobalLinSysIterative::m_coeffMatrix

protected

Definition at line 94 of file GlobalLinSysIterative.h.

Referenced by UpdateKnownSolutions().

◆ m_coeffMatrixFactor

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

protected

Definition at line 95 of file GlobalLinSysIterative.h.

Referenced by DoProjection(), and UpdateKnownSolutions().

◆ m_ipivot

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

protected

Definition at line 96 of file GlobalLinSysIterative.h.

Referenced by DoProjection(), and UpdateKnownSolutions().

◆ m_isAconjugate

bool Nektar::MultiRegions::GlobalLinSysIterative::m_isAconjugate

protected

Definition at line 98 of file GlobalLinSysIterative.h.

Referenced by GlobalLinSysIterative(), and v_SolveLinearSystem().

◆ m_linsol

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

protected

Definition at line 103 of file GlobalLinSysIterative.h.

Referenced by DoProjection(), and v_SolveLinearSystem().

◆ m_linSysIterSolver

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

protected

Iterative solver: Conjugate Gradient, GMRES.

Definition at line 89 of file GlobalLinSysIterative.h.

Referenced by GlobalLinSysIterative(), and v_SolveLinearSystem().

◆ m_map

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

protected

Global to universal unique map.

Definition at line 62 of file GlobalLinSysIterative.h.

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

◆ m_maxiter

int Nektar::MultiRegions::GlobalLinSysIterative::m_maxiter

protected

maximum iterations

Definition at line 65 of file GlobalLinSysIterative.h.

Referenced by GlobalLinSysIterative().

◆ m_NekSysOp

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

protected

Definition at line 101 of file GlobalLinSysIterative.h.

Referenced by v_SolveLinearSystem().

◆ m_numPrevSols

int Nektar::MultiRegions::GlobalLinSysIterative::m_numPrevSols

protected

Definition at line 99 of file GlobalLinSysIterative.h.

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

◆ m_numSuccessiveRHS

int Nektar::MultiRegions::GlobalLinSysIterative::m_numSuccessiveRHS

protected

Definition at line 97 of file GlobalLinSysIterative.h.

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

◆ m_precon

PreconditionerSharedPtr Nektar::MultiRegions::GlobalLinSysIterative::m_precon

protected

Definition at line 76 of file GlobalLinSysIterative.h.

◆ m_precontype

MultiRegions::PreconditionerType Nektar::MultiRegions::GlobalLinSysIterative::m_precontype

protected

Definition at line 78 of file GlobalLinSysIterative.h.

◆ m_prevBasis

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

protected

Definition at line 93 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 92 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 74 of file GlobalLinSysIterative.h.

Referenced by Set_Rhs_Magnitude(), and 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 71 of file GlobalLinSysIterative.h.

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

◆ m_root

bool Nektar::MultiRegions::GlobalLinSysIterative::m_root

protected

Root if parallel.

Definition at line 86 of file GlobalLinSysIterative.h.

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

◆ m_tolerance

NekDouble Nektar::MultiRegions::GlobalLinSysIterative::m_tolerance

protected

Tolerance of iterative solver.

Definition at line 68 of file GlobalLinSysIterative.h.

Referenced by GlobalLinSysIterative(), and v_SolveLinearSystem().

◆ m_totalIterations

int Nektar::MultiRegions::GlobalLinSysIterative::m_totalIterations

protected

Definition at line 80 of file GlobalLinSysIterative.h.

◆ m_useProjection

bool Nektar::MultiRegions::GlobalLinSysIterative::m_useProjection

protected

Whether to apply projection technique.

Definition at line 83 of file GlobalLinSysIterative.h.

Referenced by GlobalLinSysIterative(), and 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

◆ DoMatrixMultiplyFlag()

◆ DoPreconditionerFlag()

◆ DoProjection()

◆ ResetKnownSolutionsToLatestOne()

◆ Set_Rhs_Magnitude()

◆ UpdateKnownSolutions()

◆ v_DoMatrixMultiply()

◆ v_SolveLinearSystem()

◆ v_UniqueMap()

Member Data Documentation

◆ IteratSolverdef

◆ IteratSolverlookupIds

◆ m_coeffMatrix

◆ m_coeffMatrixFactor

◆ m_ipivot

◆ m_isAconjugate

◆ m_linsol

◆ m_linSysIterSolver

◆ m_map

◆ m_maxiter

◆ m_NekSysOp

◆ m_numPrevSols

◆ m_numSuccessiveRHS

◆ m_precon

◆ m_precontype

◆ m_prevBasis

◆ m_prevLinSol

◆ m_rhs_mag_sm

◆ m_rhs_magnitude

◆ m_root

◆ m_tolerance

◆ m_totalIterations

◆ m_useProjection