43 "LinSysIterSolver",
"ConjugateGradient",
46 "LinSysIterSolver",
"ConjugateGradientLoc",
67 const std::shared_ptr<AssemblyMap> &pLocToGloMap)
71 m_useProjection(false), m_numPrevSols(0)
73 m_tolerance = pLocToGloMap->GetIterativeTolerance();
78 m_expList.lock()->GetComm()->GetRowComm();
79 m_root = (vComm->GetRank()) ?
false :
true;
89 m_matrixType.find(
"AdvectionDiffusionReaction") != string::npos;
97 "Detected ConjugateGradient solver and a "
98 "Advection-Diffusion-Reaction matrix. "
99 "Switchted to a GMRES solver for this non-symmetric matrix type. "
100 "Change LinSysIterSolver to GMRES in the session file to suppress "
107 WARNINGL0(
false,
"To use A-conjugate projection, the matrix "
108 "should be symmetric positive definite.");
132 const bool isAconjugate)
134 int numIterations = 0;
139 m_linsol->SolveSystem(nGlobal, pInput, pOutput, nDir, tol);
145 m_expList.lock()->GetComm()->GetRowComm();
148 int nNonDir = nGlobal - nDir;
163 cout <<
"No iterations made"
164 <<
" using tolerance of " << tol
200 alpha.get(), n, info);
205 alpha[0] = alphaback[latest];
214 Vmath::Vcopy(nNonDir, pInput.get() + nDir, 1, pb_s.get() + nDir, 1);
225 <<
"-bases projection reduces L2-norm of RHS from "
238 numIterations =
m_linsol->SolveSystem(nGlobal, pb_s, tmpx_s, nDir, tol);
273 const bool isAconjugate)
276 int nNonDir = nGlobal - nDir;
282 m_expList.lock()->GetComm()->GetRowComm();
294 newBasis = newX + nDir;
298 newBasis = tmpAx_s + nDir;
313 tmp = tmpAx_s + nDir, 1);
317 if (i == insertLocation)
321 int skip = i > insertLocation;
335 if (i == insertLocation)
339 int iskip = i > insertLocation;
342 if (j == insertLocation)
346 int jskip = j > insertLocation;
347 tilCoeffMatrix->SetValue(i - iskip, j - jskip,
357 1, tilCoeffMatrix->GetPtr(), 1);
360 int n, info1 = 0, info2 = 0, info3 = 0;
363 n = tilCoeffMatrix->GetRows();
365 tilCoeffMatrix->GetStorageSize());
366 Vmath::Vcopy(tilCoeffMatrix->GetStorageSize(), tilCoeffMatrix->GetPtr(),
367 1, tilCoeffMatrixFactor, 1);
368 Lapack::Dsptrf(
'U', n, tilCoeffMatrixFactor.get(), ipivot.get(), info1);
372 Lapack::Dsptrs(
'U', n, 1, tilCoeffMatrixFactor.get(), ipivot.get(),
373 invMy_s.get(), n, info2);
379 y_s[0] = y_s[latest - (latest > insertLocation)];
386 for (
int i = 0; i <
m_prevBasis.size() - fullbuffer; i++)
388 residual -= y_s[i] * invMy_s[i];
392 cout <<
"SuccessiveRHS: residual " << residual;
394 if (residual < epsilon)
398 cout <<
" < " << epsilon <<
", reject" << endl;
410 1, newCoeffMatrix->GetPtr(), 1);
411 newCoeffMatrix->SetValue(insertLocation, insertLocation, 1.);
414 if (i == insertLocation)
418 int iskip = i > insertLocation;
419 newCoeffMatrix->SetValue(insertLocation, i, y_s[i - iskip]);
426 newCoeffMatrix->SetValue(insertLocation, insertLocation, 1.);
429 newCoeffMatrix->SetValue(insertLocation, i, y_s[i]);
432 newCoeffMatrix->SetValue(i, j,
m_coeffMatrix->GetValue(i, j));
436 n = newCoeffMatrix->GetRows();
438 Vmath::Vcopy(newCoeffMatrix->GetStorageSize(), newCoeffMatrix->GetPtr(), 1,
439 coeffMatrixFactor, 1);
440 Lapack::Dsptrf(
'U', n, coeffMatrixFactor.get(), ipivot.get(), info3);
445 cout <<
" >= " << epsilon <<
", reject (Dsptrf fails)" << endl;
451 cout <<
" >= " << epsilon <<
", accept" << endl;
#define WARNINGL0(condition, msg)
void DefineAssembleLoc(FuncPointerT func, ObjectPointerT obj)
void DefineNekSysLhsEval(FuncPointerT func, ObjectPointerT obj)
void DefineNekSysPrecon(FuncPointerT func, ObjectPointerT obj)
static std::string RegisterEnumValue(std::string pEnum, std::string pString, int pEnumValue)
Registers an enumeration value.
static std::string RegisterDefaultSolverInfo(const std::string &pName, const std::string &pValue)
Registers the default string value of a solver info property.
static std::shared_ptr< DataType > AllocateSharedPtr(const Args &...args)
Allocate a shared pointer from the memory pool.
const std::weak_ptr< ExpList > m_expList
Local Matrix System.
GlobalLinSysIterative(const GlobalLinSysKey &pKey, const std::weak_ptr< ExpList > &pExpList, const std::shared_ptr< AssemblyMap > &pLocToGloMap)
Constructor for full direct matrix solve.
static std::string IteratSolverdef
Array< OneD, int > m_ipivot
std::vector< Array< OneD, NekDouble > > m_prevBasis
LibUtilities::NekLinSysIterSharedPtr m_linsol
bool m_isAbsoluteTolerance
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
void UpdateKnownSolutions(const int pGlobalBndDofs, const Array< OneD, const NekDouble > &pSolution, const int pNumDirBndDofs, const bool isAconjugate)
Array< OneD, NekDouble > m_coeffMatrixFactor
void DoAssembleLocFlag(const Array< OneD, NekDouble > &pInput, Array< OneD, NekDouble > &pOutput, const bool &ZeroDir)
int ResetKnownSolutionsToLatestOne()
bool m_isNonSymmetricLinSys
void DoPreconditionerFlag(const Array< OneD, NekDouble > &pInput, Array< OneD, NekDouble > &pOutput, const bool &isLocal=false)
std::vector< Array< OneD, NekDouble > > m_prevLinSol
Storage for solutions to previous linear problems.
void DoMatrixMultiplyFlag(const Array< OneD, NekDouble > &pInput, Array< OneD, NekDouble > &pOutput, const bool &controlFlag)
static std::string IteratSolverlookupIds[]
bool m_useProjection
Whether to apply projection technique.
NekDouble m_tolerance
Tolerance of iterative solver.
Array< OneD, int > m_map
Global to universal unique map.
NekDouble m_rhs_magnitude
dot product of rhs to normalise stopping criterion
DNekMatSharedPtr m_coeffMatrix
LibUtilities::NekSysOperators m_NekSysOp
std::string m_linSysIterSolver
Iterative solver: Conjugate Gradient, GMRES.
~GlobalLinSysIterative() override
bool m_root
Root if parallel.
Describe a linear system.
StdRegions::MatrixType GetMatrixType() const
Return the matrix type.
static void Dsptrs(const char &uplo, const int &n, const int &nrhs, const double *ap, const int *ipiv, double *b, const int &ldb, int &info)
Solve a real packed-symmetric matrix problem using Bunch-Kaufman pivoting.
static void Dsptrf(const char &uplo, const int &n, double *ap, int *ipiv, int &info)
factor a real packed-symmetric matrix using Bunch-Kaufman pivoting.
std::shared_ptr< Comm > CommSharedPtr
Pointer to a Communicator object.
static PreconditionerSharedPtr NullPreconditionerSharedPtr
@ eGMRESLoc
GMRES in Local storage.
@ eConjugateGradient
Conjugate Gradient.
static const NekDouble kNekUnsetDouble
static const NekDouble kNekSparseNonZeroTol
static const NekDouble kNekZeroTol
const char *const MatrixTypeMap[]
std::shared_ptr< DNekMat > DNekMatSharedPtr
T Dot2(int n, const T *w, const T *x, const int *y)
dot product
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*x.
void Zero(int n, T *x, const int incx)
Zero vector.
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
scalarT< T > abs(scalarT< T > in)
scalarT< T > sqrt(scalarT< T > in)