42 namespace MultiRegions
51 GlobalLinSysIterative::GlobalLinSysIterative(
53 const boost::weak_ptr<ExpList> &pExpList,
54 const boost::shared_ptr<AssemblyMap>
61 m_useProjection(false),
64 m_tolerance = pLocToGloMap->GetIterativeTolerance();
65 m_maxiter = pLocToGloMap->GetMaxIterations();
68 m_root = (vComm->GetRank())?
false :
true;
72 if((successiveRHS = pLocToGloMap->GetSuccessiveRHS()))
125 =
m_expList.lock()->GetComm()->GetRowComm();
128 int nNonDir = nGlobal - nDir;
200 pInput.get() + nDir, 1,
201 pb_s.get() + nDir, 1);
232 =
m_expList.lock()->GetComm()->GetRowComm();
235 int nNonDir = nGlobal - nDir;
247 return std::sqrt(anorm_sq);
260 int nNonDir = nGlobal - nDir;
264 =
m_expList.lock()->GetComm()->GetRowComm();
292 tmp1 = newX + nDir, 1,
293 tmp2 = px_s + nDir, 1);
333 Vmath::Smul(nNonDir, 1.0/anorm, px_s.get() + nDir, 1, px_s.get() + nDir, 1);
371 =
m_expList.lock()->GetComm()->GetRowComm();
374 int nNonDir = nGlobal - nDir;
394 NekDouble alpha, beta, rho, rho_new, mu, eps, min_resid;
430 <<
" (error = " << sqrt(eps/m_rhs_magnitude)
431 <<
", rhs_mag = " << sqrt(m_rhs_magnitude) <<
")"
437 m_precon->DoPreconditioner(r_A, tmp = w_A + nDir);
471 <<
" (error = " << sqrt(eps/m_rhs_magnitude)
472 <<
", rhs_mag = " << sqrt(m_rhs_magnitude) <<
")"
476 "Exceeded maximum number of iterations");
480 Vmath::Svtvp(nNonDir, beta, &p_A[0], 1, &w_A[nDir], 1, &p_A[0], 1);
481 Vmath::Svtvp(nNonDir, beta, &q_A[0], 1, &s_A[nDir], 1, &q_A[0], 1);
484 Vmath::Svtvp(nNonDir, alpha, &p_A[0], 1, &pOutput[nDir], 1, &pOutput[nDir], 1);
487 Vmath::Svtvp(nNonDir, -alpha, &q_A[0], 1, &r_A[0], 1, &r_A[0], 1);
490 m_precon->DoPreconditioner(r_A, tmp = w_A + nDir);
515 rho_new = vExchange[0];
527 <<
" (error = " << sqrt(eps/m_rhs_magnitude)
528 <<
", rhs_mag = " << sqrt(m_rhs_magnitude) <<
")"
533 min_resid = min(min_resid, eps);
537 alpha = rho_new/(mu - rho_new*beta/alpha);
547 if (
m_map.num_elements() > 0)
550 &pIn[0],&pIn[0],&
m_map[0]);
553 m_expList.lock()->GetComm()->GetRowComm()->AllReduce(
560 NekDouble new_rhs_mag = (vExchange[0] > 1e-6)? vExchange[0] : 1.0;
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
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
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.
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
#define ROOTONLY_NEKERROR(type, msg)
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.