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)

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

Protected Member Functions inherited from Nektar::MultiRegions::GlobalLinSys
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 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)

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

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

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 63 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 94 of file GlobalLinSysIterative.cpp.

95 {

96 }

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 131 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 141 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 156 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 280 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 484 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 301 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
	)

privatepure 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
	)

privatevirtual

Solve the matrix system.

Implements Nektar::MultiRegions::GlobalLinSys.

Definition at line 101 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::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 107 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 106 of file GlobalLinSysIterative.h.

◆ m_coeffMatrix

DNekMatSharedPtr Nektar::MultiRegions::GlobalLinSysIterative::m_coeffMatrix

protected

Definition at line 95 of file GlobalLinSysIterative.h.

Referenced by UpdateKnownSolutions().

◆ m_coeffMatrixFactor

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

protected

Definition at line 96 of file GlobalLinSysIterative.h.

Referenced by DoProjection(), and UpdateKnownSolutions().

◆ m_ipivot

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

protected

Definition at line 97 of file GlobalLinSysIterative.h.

Referenced by DoProjection(), and UpdateKnownSolutions().

◆ m_isAconjugate

bool Nektar::MultiRegions::GlobalLinSysIterative::m_isAconjugate

protected

Definition at line 99 of file GlobalLinSysIterative.h.

Referenced by GlobalLinSysIterative(), and v_SolveLinearSystem().

◆ m_linsol

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

protected

Definition at line 104 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 90 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 63 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 66 of file GlobalLinSysIterative.h.

Referenced by GlobalLinSysIterative().

◆ m_NekSysOp

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

protected

Definition at line 102 of file GlobalLinSysIterative.h.

Referenced by v_SolveLinearSystem().

◆ m_numPrevSols

int Nektar::MultiRegions::GlobalLinSysIterative::m_numPrevSols

protected

Definition at line 100 of file GlobalLinSysIterative.h.

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

◆ m_numSuccessiveRHS

int Nektar::MultiRegions::GlobalLinSysIterative::m_numSuccessiveRHS

protected

Definition at line 98 of file GlobalLinSysIterative.h.

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

◆ m_precon

PreconditionerSharedPtr Nektar::MultiRegions::GlobalLinSysIterative::m_precon

protected

Definition at line 77 of file GlobalLinSysIterative.h.

◆ m_precontype

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

protected

Definition at line 79 of file GlobalLinSysIterative.h.

◆ m_prevBasis

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

protected

Definition at line 94 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 93 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 75 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 72 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 87 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 69 of file GlobalLinSysIterative.h.

Referenced by GlobalLinSysIterative(), and v_SolveLinearSystem().

◆ m_totalIterations

int Nektar::MultiRegions::GlobalLinSysIterative::m_totalIterations

protected

Definition at line 81 of file GlobalLinSysIterative.h.

◆ m_useProjection

bool Nektar::MultiRegions::GlobalLinSysIterative::m_useProjection

protected

Whether to apply projection technique.

Definition at line 84 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