Nektar::MultiRegions::GlobalLinSysIterative Class Reference

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.
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.
virtual	~GlobalLinSys ()
const GlobalLinSysKey &	GetKey (void) const
	Returns the key associated with the system.
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.
boost::shared_ptr< GlobalLinSys >	GetSharedThisPtr ()
	Returns a shared pointer to the current object.
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.

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.
void	DoConjugateGradient (const int pNumRows, const Array< OneD, const NekDouble > &pInput, Array< OneD, NekDouble > &pOutput, const AssemblyMapSharedPtr &locToGloMap, const int pNumDir)
	Actual iterative solve.
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.
virtual DNekScalMatSharedPtr	v_GetBlock (unsigned int n)
	Retrieves the block matrix from n-th expansion using the matrix key provided by the m_linSysKey.
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.
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.

Protected Attributes
Array< OneD, int >	m_map
	Global to universal unique map.
NekDouble	m_tolerance
	Tolerance of iterative solver.
NekDouble	m_rhs_magnitude
	dot product of rhs to normalise stopping criterion
PreconditionerSharedPtr	m_precon
MultiRegions::PreconditionerType	m_precontype
int	m_totalIterations
bool	m_useProjection
	Whether to apply projection technique.
bool	m_root
	Provide verbose output and root if parallel.
bool	m_verbose
boost::circular_buffer< Array < OneD, NekDouble > >	m_prevLinSol
	Storage for solutions to previous linear problems.
int	m_numPrevSols
	Total counter of previous solutions.
Protected Attributes inherited from Nektar::MultiRegions::GlobalLinSys
const GlobalLinSysKey	m_linSysKey
	Key associated with this linear system.
const boost::weak_ptr< ExpList >	m_expList
	Local Matrix System.
const map< int, RobinBCInfoSharedPtr >	m_robinBCInfo
	Robin boundary info.

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

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

                : GlobalLinSys(pKey, pExpList, pLocToGloMap),
                  m_rhs_magnitude(NekConstants::kNekUnsetDouble),
                  m_precon(NullPreconditionerSharedPtr),
                  m_totalIterations(0),
                  m_useProjection(false),
                  m_numPrevSols(0)
        {
            LibUtilities::SessionReaderSharedPtr vSession
                                            = pExpList.lock()->GetSession();

            m_tolerance = pLocToGloMap->GetIterativeTolerance();

            LibUtilities::CommSharedPtr vComm = m_expList.lock()->GetComm()->GetRowComm();
            m_root    = (vComm->GetRank())? false : true;
            m_verbose = (vSession->DefinesCmdLineArgument("verbose"))? true :false;
            
            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.

        {
        }

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::LibUtilities::NekFactory< tKey, tBase, >::CreateInstance(), Vmath::Dot2(), Nektar::eWrapper, Nektar::MultiRegions::GetPreconFactory(), Nektar::MultiRegions::GlobalLinSys::GetSharedThisPtr(), Nektar::NekConstants::kNekUnsetDouble, Nektar::MultiRegions::GlobalLinSys::m_expList, m_map, m_precon, m_rhs_magnitude, m_root, m_tolerance, m_totalIterations, m_verbose, Nektar::MultiRegions::PreconditionerTypeMap, Nektar::LibUtilities::ReduceSum, Vmath::Svtvp(), v_DoMatrixMultiply(), v_UniqueMap(), and Vmath::Zero().

Referenced by DoAconjugateProjection(), and v_SolveLinearSystem().

        {
            if (!m_precon)
            {
                MultiRegions::PreconditionerType pType
                    = plocToGloMap->GetPreconType();
                std::string PreconType
                    = MultiRegions::PreconditionerTypeMap[pType];
                v_UniqueMap();
                m_precon = GetPreconFactory().CreateInstance(
                    PreconType,GetSharedThisPtr(),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)
            {
                m_rhs_magnitude = 1.0/vExchange[2];
            }

            // 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) << ")" << endl;
                }
                m_rhs_magnitude = NekConstants::kNekUnsetDouble;
                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)
            {
                ASSERTL0(k < 5000,
                         "Exceeded maximum number of iterations (5000)");

                // 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) << ")" 
                            //<< " (bb_inv = " << sqrt(1.0/m_rhs_magnitude) << ")" 
                             << endl;
                    }
                    m_rhs_magnitude = NekConstants::kNekUnsetDouble;
                    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 534 of file GlobalLinSysIterative.cpp.

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

Referenced by 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);
            m_rhs_magnitude = (vExchange[0] > 1e-6)? vExchange[0]:1.0;
        }

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

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.

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);
            }
        }

virtual void Nektar::MultiRegions::GlobalLinSysIterative::v_UniqueMap ( )

protectedpure virtual

Referenced by DoConjugateGradient().

Member Data Documentation

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

protected

Global to universal unique map.

Definition at line 65 of file GlobalLinSysIterative.h.

Referenced by CalculateAnorm(), DoAconjugateProjection(), DoConjugateGradient(), UpdateKnownSolutions(), Nektar::MultiRegions::GlobalLinSysIterativeFull::v_UniqueMap(), and Nektar::MultiRegions::GlobalLinSysIterativeStaticCond::v_UniqueMap().

int Nektar::MultiRegions::GlobalLinSysIterative::m_numPrevSols

protected

Total counter of previous solutions.

Definition at line 90 of file GlobalLinSysIterative.h.

Referenced by DoAconjugateProjection(), and UpdateKnownSolutions().

PreconditionerSharedPtr Nektar::MultiRegions::GlobalLinSysIterative::m_precon

protected

Definition at line 73 of file GlobalLinSysIterative.h.

Referenced by DoConjugateGradient(), Nektar::MultiRegions::GlobalLinSysIterativeStaticCond::GlobalLinSysIterativeStaticCond(), Nektar::MultiRegions::GlobalLinSysIterativeStaticCond::v_AssembleSchurComplement(), Nektar::MultiRegions::GlobalLinSysIterativeStaticCond::v_BasisInvTransform(), Nektar::MultiRegions::GlobalLinSysIterativeStaticCond::v_BasisTransform(), Nektar::MultiRegions::GlobalLinSysIterativeStaticCond::v_InitObject(), and Nektar::MultiRegions::GlobalLinSysIterativeStaticCond::v_Recurse().

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

protected

Definition at line 75 of file GlobalLinSysIterative.h.

boost::circular_buffer<Array<OneD, NekDouble> > Nektar::MultiRegions::GlobalLinSysIterative::m_prevLinSol

protected

Storage for solutions to previous linear problems.

Definition at line 87 of file GlobalLinSysIterative.h.

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

NekDouble Nektar::MultiRegions::GlobalLinSysIterative::m_rhs_magnitude

protected

dot product of rhs to normalise stopping criterion

Definition at line 71 of file GlobalLinSysIterative.h.

Referenced by DoConjugateGradient(), and Set_Rhs_Magnitude().

bool Nektar::MultiRegions::GlobalLinSysIterative::m_root

protected

Provide verbose output and root if parallel.

Definition at line 83 of file GlobalLinSysIterative.h.

Referenced by DoConjugateGradient(), and GlobalLinSysIterative().

NekDouble Nektar::MultiRegions::GlobalLinSysIterative::m_tolerance

protected

Tolerance of iterative solver.

Definition at line 68 of file GlobalLinSysIterative.h.

Referenced by DoConjugateGradient(), and GlobalLinSysIterative().

int Nektar::MultiRegions::GlobalLinSysIterative::m_totalIterations

protected

Definition at line 77 of file GlobalLinSysIterative.h.

Referenced by DoConjugateGradient().

bool Nektar::MultiRegions::GlobalLinSysIterative::m_useProjection

protected

Whether to apply projection technique.

Definition at line 80 of file GlobalLinSysIterative.h.

Referenced by GlobalLinSysIterative(), and v_SolveLinearSystem().

bool Nektar::MultiRegions::GlobalLinSysIterative::m_verbose

protected

Definition at line 84 of file GlobalLinSysIterative.h.

Referenced by DoConjugateGradient(), and GlobalLinSysIterative().