Nektar::LibUtilities::NekLinSysIterCG Class Reference

#include <NekLinSysIterCG.h>

Inheritance diagram for Nektar::LibUtilities::NekLinSysIterCG:

Public Member Functions
	NekLinSysIterCG (const LibUtilities::SessionReaderSharedPtr &pSession, const LibUtilities::CommSharedPtr &vRowComm, const int nDimen, const NekSysKey &pKey)
	Constructor for full direct matrix solve. More...
	~NekLinSysIterCG () override=default
Public Member Functions inherited from Nektar::LibUtilities::NekLinSysIter
	NekLinSysIter (const LibUtilities::SessionReaderSharedPtr &pSession, const LibUtilities::CommSharedPtr &vRowComm, const int nDimen, const NekSysKey &pKey)
	~NekLinSysIter () override=default
void	SetUniversalUniqueMap (const Array< OneD, const int > &map)
void	SetNekLinSysTolerance (const NekDouble in)
void	SetNekLinSysMaxIterations (const unsigned int in)
int	GetNekLinSysTolerance ()
bool	IsLocal ()
Public Member Functions inherited from Nektar::LibUtilities::NekSys
	NekSys (const LibUtilities::SessionReaderSharedPtr &pSession, const LibUtilities::CommSharedPtr &vRowComm, const int nDimen, const NekSysKey &pKey)
virtual	~NekSys ()=default
void	InitObject ()
int	SolveSystem (const int nGlobal, const Array< OneD, const NekDouble > &pInput, Array< OneD, NekDouble > &pOutput, const int nDir=0)
const NekSysOperators &	GetSysOperators ()
void	SetSysOperators (const NekSysOperators &in)
void	SetFlagWarnings (bool in)
void	SetRhsMagnitude (const NekDouble mag)

Static Public Member Functions
static NekLinSysIterSharedPtr	create (const LibUtilities::SessionReaderSharedPtr &pSession, const LibUtilities::CommSharedPtr &vRowComm, const int nDimen, const NekSysKey &pKey)
Static Public Member Functions inherited from Nektar::LibUtilities::NekLinSysIter
static NekLinSysIterSharedPtr	CreateInstance (const LibUtilities::SessionReaderSharedPtr &pSession, const LibUtilities::CommSharedPtr &vRowComm, const int nDimen, const NekSysKey &pKey)
Static Public Member Functions inherited from Nektar::LibUtilities::NekSys
static NekSysSharedPtr	CreateInstance (const LibUtilities::SessionReaderSharedPtr &pSession, const LibUtilities::CommSharedPtr &vRowComm, const int nDimen, const NekSysKey &pKey)

Static Public Attributes
static std::string	className

Protected Member Functions
void	v_InitObject () override
int	v_SolveSystem (const int nGlobal, const Array< OneD, const NekDouble > &pInput, Array< OneD, NekDouble > &pOutput, const int nDir) override
Protected Member Functions inherited from Nektar::LibUtilities::NekLinSysIter
void	v_InitObject () override
void	SetUniversalUniqueMap ()
void	Set_Rhs_Magnitude (const Array< OneD, NekDouble > &pIn)
void	ConvergenceCheck (const Array< OneD, const NekDouble > &Residual)
Protected Member Functions inherited from Nektar::LibUtilities::NekSys
virtual void	v_InitObject ()
virtual void	v_SetSysOperators (const NekSysOperators &in)
virtual int	v_SolveSystem (const int nGlobal, const Array< OneD, const NekDouble > &pInput, Array< OneD, NekDouble > &pOutput, const int nDir)

Private Member Functions
void	DoConjugateGradient (const int pNumRows, const Array< OneD, const NekDouble > &pInput, Array< OneD, NekDouble > &pOutput, const int pNumDir)
	Actual iterative solve. More...

Friends
class	MemoryManager< NekLinSysIterCG >
	Support creation through MemoryManager. More...

Additional Inherited Members
Protected Attributes inherited from Nektar::LibUtilities::NekLinSysIter
Array< OneD, int >	m_map
	Global to universal unique map. More...
NekDouble	m_NekLinSysTolerance
int	m_NekLinSysMaxIterations
int	m_totalIterations = 0
bool	m_isLocal
Protected Attributes inherited from Nektar::LibUtilities::NekSys
LibUtilities::CommSharedPtr	m_rowComm
bool	m_converged
bool	m_root
bool	m_verbose
bool	m_FlagWarnings
int	m_SysDimen
NekSysOperators	m_operator
NekDouble	m_rhs_magnitude = NekConstants::kNekUnsetDouble

Detailed Description

Solves a linear system using iterative methods.

Definition at line 48 of file NekLinSysIterCG.h.

Constructor & Destructor Documentation

NekLinSysIterCG()

Nektar::LibUtilities::NekLinSysIterCG::NekLinSysIterCG	(	const LibUtilities::SessionReaderSharedPtr &	pSession,
		const LibUtilities::CommSharedPtr &	vRowComm,
		const int	nDimen,
		const NekSysKey &	pKey
	)

Constructor for full direct matrix solve.

Definition at line 52 of file NekLinSysIterCG.cpp.

    : NekLinSysIter(pSession, vRowComm, nDimen, pKey)
{
}

~NekLinSysIterCG()

Nektar::LibUtilities::NekLinSysIterCG::~NekLinSysIterCG ( )

overridedefault

Member Function Documentation

create()

static NekLinSysIterSharedPtr Nektar::LibUtilities::NekLinSysIterCG::create ( const LibUtilities::SessionReaderSharedPtr &  pSession, const LibUtilities::CommSharedPtr &  vRowComm, const int  nDimen, const NekSysKey &  pKey  )

inlinestatic

Definition at line 54 of file NekLinSysIterCG.h.

    {
        NekLinSysIterCGSharedPtr p =
            MemoryManager<NekLinSysIterCG>::AllocateSharedPtr(
                pSession, vRowComm, nDimen, pKey);
        p->InitObject();
        return p;
    }

References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), and CellMLToNektar.cellml_metadata::p.

DoConjugateGradient()

void Nektar::LibUtilities::NekLinSysIterCG::DoConjugateGradient ( const int  nGlobal, const Array< OneD, const NekDouble > &  pInput, Array< OneD, NekDouble > &  pOutput, const int  nDir  )

private

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 m_operator.DoNekSysLhsEval 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 91 of file NekLinSysIterCG.cpp.

{
    // 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;
 
    NekDouble alpha;
    NekDouble beta;
    NekDouble rho;
    NekDouble rho_new;
    NekDouble mu;
    NekDouble eps;
    Array<OneD, NekDouble> vExchange(3, 0.0);
 
    // Copy initial residual from input
    Vmath::Vcopy(nNonDir, pInput + nDir, 1, r_A, 1);
 
    // 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);
 
    m_rowComm->AllReduce(vExchange[2], Nektar::LibUtilities::ReduceSum);
 
    eps = vExchange[2];
 
    if (m_rhs_magnitude == NekConstants::kNekUnsetDouble)
    {
        Set_Rhs_Magnitude(pInput);
    }
 
    m_totalIterations = 0;
 
    // If input residual is less than tolerance skip solve.
    if (eps < m_NekLinSysTolerance * m_NekLinSysTolerance * m_rhs_magnitude)
    {
        if (m_verbose && m_root)
        {
            cout << "CG iterations made = " << m_totalIterations
                 << " using tolerance of " << m_NekLinSysTolerance
                 << " (error = " << sqrt(eps / m_rhs_magnitude)
                 << ", rhs_mag = " << sqrt(m_rhs_magnitude) << ")" << endl;
        }
        return;
    }
 
    m_operator.DoNekSysPrecon(r_A, tmp = w_A + nDir);
    m_operator.DoNekSysLhsEval(w_A, s_A);
 
    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);
 
    m_rowComm->AllReduce(vExchange, Nektar::LibUtilities::ReduceSum);
 
    rho               = vExchange[0];
    mu                = vExchange[1];
    beta              = 0.0;
    alpha             = rho / mu;
    m_totalIterations = 1;
 
    // Continue until convergence
    while (true)
    {
        if (m_totalIterations > m_NekLinSysMaxIterations)
        {
            if (m_root)
            {
                cout << "CG iterations made = " << m_totalIterations
                     << " using tolerance of " << m_NekLinSysTolerance
                     << " (error = " << sqrt(eps / m_rhs_magnitude)
                     << ", rhs_mag = " << sqrt(m_rhs_magnitude) << ")" << endl;
            }
            ROOTONLY_NEKERROR(ErrorUtil::efatal,
                              "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_operator.DoNekSysPrecon(r_A, tmp = w_A + nDir);
 
        // Perform the method-specific matrix-vector multiply operation.
        m_operator.DoNekSysLhsEval(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
        m_rowComm->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_NekLinSysTolerance * m_NekLinSysTolerance * m_rhs_magnitude)
        {
            if (m_verbose && m_root)
            {
                cout << "CG iterations made = " << m_totalIterations
                     << " using tolerance of " << m_NekLinSysTolerance
                     << " (error = " << sqrt(eps / m_rhs_magnitude)
                     << ", rhs_mag = " << sqrt(m_rhs_magnitude) << ")" << endl;
            }
            break;
        }
 
        // Compute search direction and solution coefficients
        beta  = rho_new / rho;
        alpha = rho_new / (mu - rho_new * beta / alpha);
        rho   = rho_new;
    }
}

References Nektar::LibUtilities::beta, Nektar::LibUtilities::NekSysOperators::DoNekSysLhsEval(), Nektar::LibUtilities::NekSysOperators::DoNekSysPrecon(), Vmath::Dot2(), Nektar::ErrorUtil::efatal, Nektar::NekConstants::kNekUnsetDouble, Nektar::LibUtilities::NekLinSysIter::m_map, Nektar::LibUtilities::NekLinSysIter::m_NekLinSysMaxIterations, Nektar::LibUtilities::NekLinSysIter::m_NekLinSysTolerance, Nektar::LibUtilities::NekSys::m_operator, Nektar::LibUtilities::NekSys::m_rhs_magnitude, Nektar::LibUtilities::NekSys::m_root, Nektar::LibUtilities::NekSys::m_rowComm, Nektar::LibUtilities::NekLinSysIter::m_totalIterations, Nektar::LibUtilities::NekSys::m_verbose, Nektar::LibUtilities::ReduceSum, ROOTONLY_NEKERROR, Nektar::LibUtilities::NekLinSysIter::Set_Rhs_Magnitude(), tinysimd::sqrt(), Vmath::Svtvp(), Vmath::Vcopy(), and Vmath::Zero().

Referenced by v_SolveSystem().

v_InitObject()

void Nektar::LibUtilities::NekLinSysIterCG::v_InitObject ( )

overrideprotectedvirtual

Reimplemented from Nektar::LibUtilities::NekLinSysIter.

Definition at line 60 of file NekLinSysIterCG.cpp.

{
    NekLinSysIter::v_InitObject();
}

References Nektar::LibUtilities::NekLinSysIter::v_InitObject().

v_SolveSystem()

int Nektar::LibUtilities::NekLinSysIterCG::v_SolveSystem ( const int  nGlobal, const Array< OneD, const NekDouble > &  pInput, Array< OneD, NekDouble > &  pOutput, const int  nDir  )

overrideprotectedvirtual

Reimplemented from Nektar::LibUtilities::NekSys.

Definition at line 68 of file NekLinSysIterCG.cpp.

{
    DoConjugateGradient(nGlobal, pInput, pOutput, nDir);
 
    return m_totalIterations;
}

References DoConjugateGradient(), and Nektar::LibUtilities::NekLinSysIter::m_totalIterations.

Friends And Related Function Documentation

MemoryManager< NekLinSysIterCG >

friend class MemoryManager< NekLinSysIterCG >

friend

Support creation through MemoryManager.

Definition at line 46 of file NekLinSysIterCG.h.

Member Data Documentation

className

string Nektar::LibUtilities::NekLinSysIterCG::className

static

Initial value:

=
    LibUtilities::GetNekLinSysIterFactory().RegisterCreatorFunction(
        "ConjugateGradient", NekLinSysIterCG::create,
        "NekLinSysIterCG solver.")

Definition at line 66 of file NekLinSysIterCG.h.