42 namespace MultiRegions
51 GlobalLinSysIterative::GlobalLinSysIterative(
53 const boost::weak_ptr<ExpList> &pExpList,
54 const boost::shared_ptr<AssemblyMap>
61 m_useProjection(false),
65 = pExpList.lock()->GetSession();
67 m_tolerance = pLocToGloMap->GetIterativeTolerance();
68 m_maxiter = pLocToGloMap->GetMaxIterations();
71 m_root = (vComm->GetRank())?
false :
true;
72 m_verbose = (vSession->DefinesCmdLineArgument(
"verbose"))?
true :
false;
76 if((successiveRHS = pLocToGloMap->GetSuccessiveRHS()))
129 =
m_expList.lock()->GetComm()->GetRowComm();
132 int nNonDir = nGlobal - nDir;
204 pInput.get() + nDir, 1,
205 pb_s.get() + nDir, 1);
236 =
m_expList.lock()->GetComm()->GetRowComm();
239 int nNonDir = nGlobal - nDir;
251 return std::sqrt(anorm_sq);
264 int nNonDir = nGlobal - nDir;
268 =
m_expList.lock()->GetComm()->GetRowComm();
296 tmp1 = newX + nDir, 1,
297 tmp2 = px_s + nDir, 1);
337 Vmath::Smul(nNonDir, 1.0/anorm, px_s.get() + nDir, 1, px_s.get() + nDir, 1);
375 =
m_expList.lock()->GetComm()->GetRowComm();
378 int nNonDir = nGlobal - nDir;
398 NekDouble alpha, beta, rho, rho_new, mu, eps, min_resid;
432 <<
" (error = " << sqrt(eps/m_rhs_magnitude)
433 <<
", rhs_mag = " << sqrt(m_rhs_magnitude) <<
")"
440 m_precon->DoPreconditioner(r_A, tmp = w_A + nDir);
473 <<
" (error = " << sqrt(eps/m_rhs_magnitude)
474 <<
", rhs_mag = " << sqrt(m_rhs_magnitude) <<
")"
477 "Exceeded maximum number of iterations");
482 Vmath::Svtvp(nNonDir, beta, &p_A[0], 1, &w_A[nDir], 1, &p_A[0], 1);
483 Vmath::Svtvp(nNonDir, beta, &q_A[0], 1, &s_A[nDir], 1, &q_A[0], 1);
486 Vmath::Svtvp(nNonDir, alpha, &p_A[0], 1, &pOutput[nDir], 1, &pOutput[nDir], 1);
489 Vmath::Svtvp(nNonDir, -alpha, &q_A[0], 1, &r_A[0], 1, &r_A[0], 1);
492 m_precon->DoPreconditioner(r_A, tmp = w_A + nDir);
517 rho_new = vExchange[0];
529 <<
" (error = " << sqrt(eps/m_rhs_magnitude)
530 <<
", rhs_mag = " << sqrt(m_rhs_magnitude) <<
")"
535 min_resid = min(min_resid, eps);
539 alpha = rho_new/(mu - rho_new*beta/alpha);
551 m_expList.lock()->GetComm()->GetRowComm()->AllReduce(
558 NekDouble new_rhs_mag = (vExchange[0] > 1e-6)? vExchange[0] : 1.0;
#define ASSERTL0(condition, msg)
boost::shared_ptr< AssemblyMap > AssemblyMapSharedPtr
NekDouble m_rhs_mag_sm
cnt to how many times rhs_magnitude is called
void Set_Rhs_Magnitude(const NekVector< NekDouble > &pIn)
void UpdateKnownSolutions(const int pGlobalBndDofs, const Array< OneD, const NekDouble > &pSolution, const int pNumDirBndDofs)
bool m_root
Provide verbose output and root if parallel.
void Svtvp(int n, const T alpha, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
svtvp (scalar times vector plus vector): z = alpha*x + y
bool m_useProjection
Whether to apply projection technique.
boost::circular_buffer< Array< OneD, NekDouble > > m_prevLinSol
Storage for solutions to previous linear problems.
virtual void v_DoMatrixMultiply(const Array< OneD, NekDouble > &pInput, Array< OneD, NekDouble > &pOutput)=0
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
NekDouble CalculateAnorm(const int nGlobal, const Array< OneD, const NekDouble > &in, const int nDir)
PreconditionerSharedPtr m_precon
boost::shared_ptr< Comm > CommSharedPtr
Pointer to a Communicator object.
static const NekDouble kNekZeroTol
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*y.
static PreconditionerSharedPtr NullPreconditionerSharedPtr
NekDouble m_tolerance
Tolerance of iterative solver.
int m_numPrevSols
Total counter of previous solutions.
Array< OneD, int > m_map
Global to universal unique map.
static const NekDouble kNekUnsetDouble
Describe a linear system.
T Dot(int n, const T *w, const T *x)
vvtvp (vector times vector times vector): z = w*x*y
unsigned int GetDimension() const
Returns the number of dimensions for the point.
void DoAconjugateProjection(const int pNumRows, const Array< OneD, const NekDouble > &pInput, Array< OneD, NekDouble > &pOutput, const AssemblyMapSharedPtr &locToGloMap, const int pNumDir)
A-conjugate projection technique.
NekDouble m_rhs_magnitude
dot product of rhs to normalise stopping criterion
int m_maxiter
maximum iterations
T Dot2(int n, const T *w, const T *x, const int *y)
vvtvp (vector times vector times vector): z = w*x*y
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.
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)
void DoConjugateGradient(const int pNumRows, const Array< OneD, const NekDouble > &pInput, Array< OneD, NekDouble > &pOutput, const AssemblyMapSharedPtr &locToGloMap, const int pNumDir)
Actual iterative solve.
PreconditionerSharedPtr CreatePrecon(AssemblyMapSharedPtr asmMap)
Create a preconditioner object from the parameters defined in the supplied assembly map...
virtual ~GlobalLinSysIterative()
virtual void v_UniqueMap()=0
const boost::weak_ptr< ExpList > m_expList
Local Matrix System.