#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	setRhsMagnitude (const NekDouble mag)

void	SetNekLinSysMaxIterations (const unsigned int in)

void	SetLinSysMaxStorage (const unsigned int in)

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 ()

void	SetSysOperators (const NekSysOperators &in)

const NekSysOperators &	GetSysOperators ()

int	SolveSystem (const int nGlobal, const Array< OneD, const NekDouble > &pInput, Array< OneD, NekDouble > &pOutput, const int nDir, const NekDouble tol=1.0E-7, const NekDouble factor=1.0)

bool	ConvergenceCheck (const int nIteration, const Array< OneD, const NekDouble > &Residual, const NekDouble tol=1.0E-7)

void	SetFlagWarnings (bool in)

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, const NekDouble tol, const NekDouble factor) override

Protected Member Functions inherited from Nektar::LibUtilities::NekLinSysIter
void	Set_Rhs_Magnitude (const NekVector< NekDouble > &pIn)

void	Set_Rhs_Magnitude (const Array< OneD, NekDouble > &pIn)

void	SetUniversalUniqueMap ()

void	v_InitObject () override

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, const NekDouble tol, const NekDouble factor)

virtual bool	v_ConvergenceCheck (const int nIteration, const Array< OneD, const NekDouble > &Residual, const NekDouble tol)

virtual void	v_NekSysInitialGuess (const Array< OneD, const NekDouble > &pInput, Array< OneD, NekDouble > &pguess)

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_rhs_magnitude = NekConstants::kNekUnsetDouble
	Dot product of rhs to normalise stopping criterion. More...

int	m_totalIterations = 0

NekDouble	m_prec_factor = 1.0

int	m_LinSysMaxStorage

bool	m_isLocal

Protected Attributes inherited from Nektar::LibUtilities::NekSys
int	m_maxiter
	Maximum iterations. More...

NekDouble	m_tolerance
	Tolerance of iterative solver. More...

LibUtilities::CommSharedPtr	m_rowComm
	Communicate. More...

bool	m_converged
	Whether the iteration has been converged. More...

bool	m_root
	Root if parallel. More...

bool	m_verbose
	Verbose. More...

bool	m_FlagWarnings

NekSysOperators	m_operator
	Operators. More...

int	m_SysDimen
	The dimension of the system. More...

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 94 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)
    {
        NekVector<NekDouble> inGlob(nGlobal, pInput, eWrapper);
        Set_Rhs_Magnitude(inGlob);
    }
 
    m_totalIterations = 0;
 
    // 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_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_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;
            }
            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_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;
        }
 
        // Compute search direction and solution coefficients
        beta  = rho_new / rho;
        alpha = rho_new / (mu - rho_new * beta / alpha);
        rho   = rho_new;
    }
}

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,
		const NekDouble	tol,
		const NekDouble	factor
	)

overrideprotectedvirtual

Reimplemented from Nektar::LibUtilities::NekSys.

Definition at line 68 of file NekLinSysIterCG.cpp.

{
    m_tolerance = max(tol, 1.0E-16);
 
    DoConjugateGradient(nGlobal, pInput, pOutput, nDir);
 
    return m_totalIterations;
}

References DoConjugateGradient(), Nektar::LibUtilities::NekSys::m_tolerance, 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 65 of file NekLinSysIterCG.h.

Public Member Functions

Static Public Member Functions

Static Public Attributes

Protected Member Functions

Private Member Functions

Friends

Additional Inherited Members