A global linear system. More...

#include <GlobalLinSysIterative.h>

Inheritance diagram for Nektar::MultiRegions::GlobalLinSysIterative:

Collaboration diagram for Nektar::MultiRegions::GlobalLinSysIterative:

Public Member Functions
	GlobalLinSysIterative (const GlobalLinSysKey &pKey, const boost::weak_ptr< ExpList > &pExpList, const boost::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 boost::weak_ptr< ExpList > &pExpList, const boost::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 boost::weak_ptr< ExpList > &	GetLocMat (void) const

void	InitObject ()

void	Initialise (const boost::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...

boost::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	DoAconjugateProjection (const int pNumRows, const Array< OneD, const NekDouble > &pInput, Array< OneD, NekDouble > &pOutput, const AssemblyMapSharedPtr &locToGloMap, const int pNumDir)
	A-conjugate projection technique. More...

void	DoConjugateGradient (const int pNumRows, const Array< OneD, const NekDouble > &pInput, Array< OneD, NekDouble > &pOutput, const AssemblyMapSharedPtr &locToGloMap, const int pNumDir)
	Actual iterative solve. 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...

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

int	m_numPrevSols
	Total counter of previous solutions. More...

Protected Attributes inherited from Nektar::MultiRegions::GlobalLinSys
const GlobalLinSysKey	m_linSysKey
	Key associated with this linear system. More...

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

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

bool	m_verbose

Private Member Functions
void	UpdateKnownSolutions (const int pGlobalBndDofs, const Array< OneD, const NekDouble > &pSolution, const int pNumDirBndDofs)

NekDouble	CalculateAnorm (const int nGlobal, const Array< OneD, const NekDouble > &in, const int nDir)

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 52 of file GlobalLinSysIterative.h.

Constructor & Destructor Documentation

Nektar::MultiRegions::GlobalLinSysIterative::GlobalLinSysIterative	(	const GlobalLinSysKey &	pKey,
		const boost::weak_ptr< ExpList > &	pExpList,
		const boost::shared_ptr< AssemblyMap > &	pLocToGloMap
	)

Constructor for full direct matrix solve.

Definition at line 51 of file GlobalLinSysIterative.cpp.

References Nektar::MultiRegions::GlobalLinSys::m_expList, m_maxiter, m_prevLinSol, m_root, m_tolerance, and m_useProjection.

                 : 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();
 
             LibUtilities::CommSharedPtr vComm = m_expList.lock()->GetComm()->GetRowComm();
             m_root    = (vComm->GetRank())? false : true;
 
             int successiveRHS;
             
             if((successiveRHS = pLocToGloMap->GetSuccessiveRHS()))
             {
                 m_prevLinSol.set_capacity(successiveRHS);
                 m_useProjection = true;
             }
             else
             {
                 m_useProjection = false;
             }
         }

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

virtual

Definition at line 83 of file GlobalLinSysIterative.cpp.

84 {

85 }

Member Function Documentation

NekDouble Nektar::MultiRegions::GlobalLinSysIterative::CalculateAnorm	(	const int	nGlobal,
		const Array< OneD, const NekDouble > &	in,
		const int	nDir
	)

private

Calculating A-norm of an input vector, A-norm(x) := sqrt( < x, Ax > )

Definition at line 225 of file GlobalLinSysIterative.cpp.

References Vmath::Dot2(), Nektar::MultiRegions::GlobalLinSys::m_expList, m_map, Nektar::LibUtilities::ReduceSum, and v_DoMatrixMultiply().

Referenced by UpdateKnownSolutions().

         {
             // Get the communicator for performing data exchanges
             LibUtilities::CommSharedPtr vComm
                 = m_expList.lock()->GetComm()->GetRowComm();
 
             // Get vector sizes
             int nNonDir = nGlobal - nDir;
 
             // Allocate array storage
             Array<OneD, NekDouble> tmpAx_s    (nGlobal, 0.0);
 
             v_DoMatrixMultiply(in, tmpAx_s);
 
             NekDouble anorm_sq = Vmath::Dot2(nNonDir,
                                              in      + nDir,
                                              tmpAx_s + nDir,
                                              m_map   + nDir);
             vComm->AllReduce(anorm_sq, Nektar::LibUtilities::ReduceSum);
             return std::sqrt(anorm_sq);
         }

void Nektar::MultiRegions::GlobalLinSysIterative::DoAconjugateProjection	(	const int	nGlobal,
		const Array< OneD, const NekDouble > &	pInput,
		Array< OneD, NekDouble > &	pOutput,
		const AssemblyMapSharedPtr &	plocToGloMap,
		const int	nDir
	)

protected

A-conjugate projection technique.

This method implements A-conjugate projection technique 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 116 of file GlobalLinSysIterative.cpp.

References DoConjugateGradient(), Vmath::Dot2(), Nektar::eWrapper, Nektar::NekConstants::kNekZeroTol, Nektar::MultiRegions::GlobalLinSys::m_expList, m_map, m_numPrevSols, m_prevLinSol, Nektar::LibUtilities::ReduceSum, UpdateKnownSolutions(), v_DoMatrixMultiply(), Vmath::Vcopy(), and Vmath::Zero().

Referenced by v_SolveLinearSystem().

         {
             // Get the communicator for performing data exchanges
             LibUtilities::CommSharedPtr vComm
                                 = m_expList.lock()->GetComm()->GetRowComm();
 
             // Get vector sizes
             int nNonDir = nGlobal - nDir;
             Array<OneD, NekDouble> tmp;
 
             if (0 == m_numPrevSols)
             {
                 // no previous solutions found, call CG
 
                 DoConjugateGradient(nGlobal, pInput, pOutput, plocToGloMap, nDir);
 
                 UpdateKnownSolutions(nGlobal, pOutput, nDir);
             }
             else
             {
                 // Create NekVector wrappers for linear algebra operations
                 NekVector<NekDouble> b     (nNonDir, pInput  + nDir, eWrapper);
                 NekVector<NekDouble> x     (nNonDir, tmp = pOutput + nDir, eWrapper);
 
                 // check the input vector (rhs) is not zero
 
                 NekDouble rhsNorm = Vmath::Dot2(nNonDir,
                                                 pInput + nDir,
                                                 pInput + nDir,
                                                 m_map + nDir);
 
                 vComm->AllReduce(rhsNorm, Nektar::LibUtilities::ReduceSum);
 
                 if (rhsNorm < NekConstants::kNekZeroTol)
                 {
                     Array<OneD, NekDouble> tmp = pOutput+nDir;
                     Vmath::Zero(nNonDir, tmp, 1);
                     return;
                 }
 
                 // 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_prevLinSol.size(), 0.0);
                 for (int i = 0; i < m_prevLinSol.size(); i++)
                 {
                     alpha[i] = Vmath::Dot2(nNonDir,
                                            m_prevLinSol[i],
                                            pInput + nDir,
                                            m_map + nDir);
                 }
                 vComm->AllReduce(alpha, Nektar::LibUtilities::ReduceSum);
 
                 for (int i = 0; i < m_prevLinSol.size(); i++)
                 {
                     if (alpha[i] < NekConstants::kNekZeroTol)
                     {
                         continue;
                     }
 
                     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);
 
                 v_DoMatrixMultiply(px_s, tmpAx_s);
 
                 pb -= tmpAx;
 
 
                 // solve the system with projected rhs
                 DoConjugateGradient(nGlobal, pb_s, tmpx_s, plocToGloMap, nDir);
 
 
                 // remainder solution + projection of previous solutions
                 x = tmpx + px;
 
                 // save the auxiliary solution to prev. known solutions
                 UpdateKnownSolutions(nGlobal, tmpx_s, nDir);
             }
         }

void Nektar::MultiRegions::GlobalLinSysIterative::DoConjugateGradient	(	const int	nGlobal,
		const Array< OneD, const NekDouble > &	pInput,
		Array< OneD, NekDouble > &	pOutput,
		const AssemblyMapSharedPtr &	plocToGloMap,
		const int	nDir
	)

protected

Actual iterative solve.

Solve a global linear system using the conjugate gradient method. We solve only for the non-Dirichlet modes. The operator is evaluated using an auxiliary function v_DoMatrixMultiply defined by the specific solver. Distributed math routines are used to support parallel execution of the solver. The implemented algorithm uses a reduced-communication reordering of the standard PCG method (Demmel, Heath and Vorst, 1993)

Parameters

pInput	Input residual of all DOFs.
pOutput	Solution vector of all DOFs.

Definition at line 355 of file GlobalLinSysIterative.cpp.

References ASSERTL0, Nektar::MultiRegions::GlobalLinSys::CreatePrecon(), Vmath::Dot2(), Nektar::eWrapper, Nektar::NekConstants::kNekUnsetDouble, Nektar::MultiRegions::GlobalLinSys::m_expList, m_map, m_maxiter, m_precon, m_rhs_magnitude, m_root, m_tolerance, m_totalIterations, Nektar::MultiRegions::GlobalLinSys::m_verbose, Nektar::LibUtilities::ReduceSum, Set_Rhs_Magnitude(), Vmath::Svtvp(), v_DoMatrixMultiply(), v_UniqueMap(), and Vmath::Zero().

Referenced by DoAconjugateProjection(), and v_SolveLinearSystem().

         {
             if (!m_precon)
             {
                 v_UniqueMap();
                 m_precon = CreatePrecon(plocToGloMap);
                 m_precon->BuildPreconditioner();
             }
 
             // Get the communicator for performing data exchanges
             LibUtilities::CommSharedPtr vComm
                 = m_expList.lock()->GetComm()->GetRowComm();
 
             // Get vector sizes
             int nNonDir = nGlobal - nDir;
 
             // Allocate array storage
             Array<OneD, NekDouble> w_A    (nGlobal, 0.0);
             Array<OneD, NekDouble> s_A    (nGlobal, 0.0);
             Array<OneD, NekDouble> p_A    (nNonDir, 0.0);
             Array<OneD, NekDouble> r_A    (nNonDir, 0.0);
             Array<OneD, NekDouble> q_A    (nNonDir, 0.0);
             Array<OneD, NekDouble> tmp;
 
             // Create NekVector wrappers for linear algebra operations
             NekVector<NekDouble> in (nNonDir,pInput  + nDir,      eWrapper);
             NekVector<NekDouble> out(nNonDir,tmp = pOutput + nDir,eWrapper);
             NekVector<NekDouble> w  (nNonDir,tmp = w_A + nDir,    eWrapper);
             NekVector<NekDouble> s  (nNonDir,tmp = s_A + nDir,    eWrapper);
             NekVector<NekDouble> p  (nNonDir,p_A,                 eWrapper);
             NekVector<NekDouble> r  (nNonDir,r_A,                 eWrapper);
             NekVector<NekDouble> q  (nNonDir,q_A,                 eWrapper);
 
             int k;
             NekDouble alpha, beta, rho, rho_new, mu, eps,  min_resid;
             Array<OneD, NekDouble> vExchange(3,0.0);
 
             // Copy initial residual from input
             r = in;
             // zero homogeneous out array ready for solution updates
             // Should not be earlier in case input vector is same as
             // output and above copy has been peformed
             Vmath::Zero(nNonDir,tmp = pOutput + nDir,1);
 
 
             // evaluate initial residual error for exit check
             vExchange[2] = Vmath::Dot2(nNonDir,
                                        r_A,
                                        r_A,
                                        m_map + nDir);
 
             vComm->AllReduce(vExchange, Nektar::LibUtilities::ReduceSum);
             
             eps       = vExchange[2];
 
             if(m_rhs_magnitude == NekConstants::kNekUnsetDouble)
             {
                 NekVector<NekDouble> inGlob (nGlobal, pInput, eWrapper);
                 Set_Rhs_Magnitude(inGlob);
             }
 
             // If input residual is less than tolerance skip solve.
             if (eps < m_tolerance * m_tolerance * m_rhs_magnitude)
             {
                 if (m_verbose && m_root)
                 {
                     cout << "CG iterations made = " << m_totalIterations 
                          << " using tolerance of "  << m_tolerance 
                          << " (error = " << sqrt(eps/m_rhs_magnitude) 
                          << ", rhs_mag = " << sqrt(m_rhs_magnitude) <<  ")" 
                          << endl;
                 }
                 return;
             }
 
             m_totalIterations = 1;
             m_precon->DoPreconditioner(r_A, tmp = w_A + nDir);
 
             v_DoMatrixMultiply(w_A, s_A);
 
             k = 0;
 
             vExchange[0] = Vmath::Dot2(nNonDir,
                                        r_A,
                                        w_A + nDir,
                                        m_map + nDir);
 
             vExchange[1] = Vmath::Dot2(nNonDir,
                                        s_A + nDir,
                                        w_A + nDir,
                                        m_map + nDir);
 
             vComm->AllReduce(vExchange, Nektar::LibUtilities::ReduceSum);
 
             rho       = vExchange[0];
             mu        = vExchange[1];
             min_resid = m_rhs_magnitude;
             beta      = 0.0;
             alpha     = rho/mu;
 
             // Continue until convergence
             while (true)
             {
                 if(k >= m_maxiter)
                 {
                     if (m_root)
                     {
                         cout << "CG iterations made = " << m_totalIterations 
                              << " using tolerance of "  << m_tolerance 
                              << " (error = " << sqrt(eps/m_rhs_magnitude)
                              << ", rhs_mag = " << sqrt(m_rhs_magnitude) <<  ")"
                              << endl;
                     }
                     ASSERTL0(false,
                              "Exceeded maximum number of iterations");
                 }
 
                 // Compute new search direction p_k, q_k
                 Vmath::Svtvp(nNonDir, beta, &p_A[0], 1, &w_A[nDir], 1, &p_A[0], 1);
                 Vmath::Svtvp(nNonDir, beta, &q_A[0], 1, &s_A[nDir], 1, &q_A[0], 1);
 
                 // Update solution x_{k+1}
                 Vmath::Svtvp(nNonDir, alpha, &p_A[0], 1, &pOutput[nDir], 1, &pOutput[nDir], 1);
 
                 // Update residual vector r_{k+1}
                 Vmath::Svtvp(nNonDir, -alpha, &q_A[0], 1, &r_A[0], 1, &r_A[0], 1);
 
                 // Apply preconditioner
                 m_precon->DoPreconditioner(r_A, tmp = w_A + nDir);
 
                 // Perform the method-specific matrix-vector multiply operation.
                 v_DoMatrixMultiply(w_A, s_A);
 
                 // <r_{k+1}, w_{k+1}>
                 vExchange[0] = Vmath::Dot2(nNonDir,
                                            r_A,
                                            w_A + nDir,
                                            m_map + nDir);
                 // <s_{k+1}, w_{k+1}>
                 vExchange[1] = Vmath::Dot2(nNonDir,
                                            s_A + nDir,
                                            w_A + nDir,
                                            m_map + nDir);
 
                 // <r_{k+1}, r_{k+1}>
                 vExchange[2] = Vmath::Dot2(nNonDir,
                                            r_A,
                                            r_A,
                                            m_map + nDir);
 
                 // Perform inner-product exchanges
                 vComm->AllReduce(vExchange, Nektar::LibUtilities::ReduceSum);
 
                 rho_new = vExchange[0];
                 mu      = vExchange[1];
                 eps     = vExchange[2];
 
                 m_totalIterations++;
                 // test if norm is within tolerance
                 if (eps < m_tolerance * m_tolerance * m_rhs_magnitude)
                 {
                     if (m_verbose && m_root)
                     {
                         cout << "CG iterations made = " << m_totalIterations 
                              << " using tolerance of "  << m_tolerance 
                              << " (error = " << sqrt(eps/m_rhs_magnitude)
                              << ", rhs_mag = " << sqrt(m_rhs_magnitude) <<  ")"
                              << endl;
                     }
                     break;
                 }
                 min_resid = min(min_resid, eps);
 
                 // Compute search direction and solution coefficients
                 beta  = rho_new/rho;
                 alpha = rho_new/(mu - rho_new*beta/alpha);
                 rho   = rho_new;
                 k++;
             }
         }

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

protected

Definition at line 541 of file GlobalLinSysIterative.cpp.

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

Referenced by DoConjugateGradient(), and Nektar::MultiRegions::GlobalLinSysIterativeStaticCond::v_PreSolve().

         {
             Array<OneD, NekDouble> vExchange(1);
             vExchange[0] = Vmath::Dot(pIn.GetDimension(),&pIn[0],&pIn[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); 
             }
         }

void Nektar::MultiRegions::GlobalLinSysIterative::UpdateKnownSolutions	(	const int	nGlobal,
		const Array< OneD, const NekDouble > &	newX,
		const int	nDir
	)

private

Updates the storage of previously known solutions. Performs normalisation of input vector wrt A-norm.

Definition at line 254 of file GlobalLinSysIterative.cpp.

References CalculateAnorm(), Vmath::Dot2(), Nektar::eWrapper, Nektar::NekConstants::kNekZeroTol, Nektar::MultiRegions::GlobalLinSys::m_expList, m_map, m_numPrevSols, m_prevLinSol, Nektar::LibUtilities::ReduceSum, Vmath::Smul(), v_DoMatrixMultiply(), and Vmath::Vcopy().

Referenced by DoAconjugateProjection().

         {
             // Get vector sizes
             int nNonDir = nGlobal - nDir;
 
             // Get the communicator for performing data exchanges
             LibUtilities::CommSharedPtr vComm
                 = m_expList.lock()->GetComm()->GetRowComm();
 
             // Check the solution is non-zero
             NekDouble solNorm = Vmath::Dot2(nNonDir,
                                             newX + nDir,
                                             newX + nDir,
                                             m_map + nDir);
             vComm->AllReduce(solNorm, Nektar::LibUtilities::ReduceSum);
 
             if (solNorm < NekConstants::kNekZeroTol)
             {
                 return;
             }
 
 
             // Allocate array storage
             Array<OneD, NekDouble> tmpAx_s    (nGlobal, 0.0);
             Array<OneD, NekDouble> px_s       (nGlobal, 0.0);
             Array<OneD, NekDouble> tmp1, tmp2;
 
             // Create NekVector wrappers for linear algebra operations
             NekVector<NekDouble> px           (nNonDir, tmp1 = px_s    + nDir, eWrapper);
             NekVector<NekDouble> tmpAx        (nNonDir, tmp2 = tmpAx_s + nDir, eWrapper);
 
 
             // calculating \tilda{x} - sum \alpha_i\tilda{x}_i
 
             Vmath::Vcopy(nNonDir,
                          tmp1 = newX + nDir, 1,
                          tmp2 = px_s + nDir, 1);
 
             if (m_prevLinSol.size() > 0)
             {
                 v_DoMatrixMultiply(newX, tmpAx_s);
             }
 
             Array<OneD, NekDouble> alpha (m_prevLinSol.size(), 0.0);
             for (int i = 0; i < m_prevLinSol.size(); i++)
             {
                 alpha[i] = Vmath::Dot2(nNonDir,
                                        m_prevLinSol[i],
                                        tmpAx_s + nDir,
                                        m_map + nDir);
             }
             vComm->AllReduce(alpha, Nektar::LibUtilities::ReduceSum);
 
             for (int i = 0; i < m_prevLinSol.size(); i++)
             {
                 if (alpha[i] < NekConstants::kNekZeroTol)
                 {
                     continue;
                 }
 
                 NekVector<NekDouble> xi (nNonDir, m_prevLinSol[i], eWrapper);
                 px -= alpha[i] * xi;
             }
 
 
             // Some solutions generated by CG are identical zeros, see
             // solutions generated for Test_Tet_equitri.xml (IncNavierStokesSolver).
             // Not going to store identically zero solutions.
 
             NekDouble anorm = CalculateAnorm(nGlobal, px_s, nDir);
             if (anorm < NekConstants::kNekZeroTol)
             {
                 return;
             }
 
             // normalisation of new solution
             Vmath::Smul(nNonDir, 1.0/anorm, px_s.get() + nDir, 1, px_s.get() + nDir, 1);
 
             // updating storage with non-Dirichlet-dof part of new solution vector
             m_prevLinSol.push_back(px_s + nDir);
             m_numPrevSols++;
         }

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 CalculateAnorm(), DoAconjugateProjection(), DoConjugateGradient(), and UpdateKnownSolutions().

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

References DoAconjugateProjection(), DoConjugateGradient(), and m_useProjection.

         {
             if (m_useProjection)
             {
                 DoAconjugateProjection(nGlobal, pInput, pOutput, plocToGloMap, nDir);
             }
             else
             {
                 // applying plain Conjugate Gradient
                 DoConjugateGradient(nGlobal, pInput, pOutput, plocToGloMap, nDir);
             }
         }