#include <NekLinSysIterGMRES.h>

Inheritance diagram for Nektar::LibUtilities::NekLinSysIterGMRES:

Public Member Functions
	NekLinSysIterGMRES (const LibUtilities::SessionReaderSharedPtr &pSession, const LibUtilities::CommSharedPtr &vRowComm, const int nDimen, const NekSysKey &pKey=NekSysKey())

	~NekLinSysIterGMRES () override=default

int	GetMaxLinIte ()

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)

Protected Attributes
int	m_maxrestart

int	m_KrylovMaxHessMatBand

bool	m_NekLinSysLeftPrecon = false

bool	m_NekLinSysRightPrecon = true

bool	m_DifferenceFlag0 = false

bool	m_DifferenceFlag1 = false

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...

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

NekDouble	DoGmresRestart (const bool restarted, const bool truncted, const int nGlobal, const Array< OneD, const NekDouble > &pInput, Array< OneD, NekDouble > &pOutput, const int nDir)
	Actual iterative gmres solver for one restart. More...

void	DoArnoldi (const int starttem, const int endtem, const int nGlobal, const int nDir, Array< OneD, NekDouble > &w, Array< OneD, NekDouble > &Vsingle1, Array< OneD, NekDouble > &Vsingle2, Array< OneD, NekDouble > &hsingle)

void	DoGivensRotation (const int starttem, const int endtem, const int nGlobal, const int nDir, Array< OneD, NekDouble > &c, Array< OneD, NekDouble > &s, Array< OneD, NekDouble > &hsingle, Array< OneD, NekDouble > &eta)

void	DoBackward (const int number, Array< OneD, Array< OneD, NekDouble > > &A, const Array< OneD, const NekDouble > &b, Array< OneD, NekDouble > &y)

Private Attributes
Array< OneD, Array< OneD, NekDouble > >	m_hes

Array< OneD, Array< OneD, NekDouble > >	m_Upper

Array< OneD, Array< OneD, NekDouble > >	m_V_total

Static Private Attributes
static std::string	lookupIds []

static std::string	def

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

Detailed Description

Solves a linear system using iterative methods.

Definition at line 46 of file NekLinSysIterGMRES.h.

Constructor & Destructor Documentation

◆ NekLinSysIterGMRES()

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

Definition at line 52 of file NekLinSysIterGMRES.cpp.

    : NekLinSysIter(pSession, vRowComm, nDimen, pKey)
{
    m_NekLinSysLeftPrecon  = pKey.m_NekLinSysLeftPrecon;
    m_NekLinSysRightPrecon = pKey.m_NekLinSysRightPrecon;
 
    m_KrylovMaxHessMatBand = pKey.m_KrylovMaxHessMatBand;
 
    m_maxrestart = ceil(NekDouble(m_maxiter) / NekDouble(m_LinSysMaxStorage));
    m_LinSysMaxStorage = min(m_maxiter, m_LinSysMaxStorage);
 
    m_DifferenceFlag0 = pKey.m_DifferenceFlag0;
    m_DifferenceFlag1 = pKey.m_DifferenceFlag1;
 
    // Allocate array storage of coefficients
    // Hessenburg matrix
    m_hes = Array<OneD, Array<OneD, NekDouble>>(m_LinSysMaxStorage);
    for (size_t nd = 0; nd < m_LinSysMaxStorage; nd++)
    {
        m_hes[nd] = Array<OneD, NekDouble>(m_LinSysMaxStorage + 1, 0.0);
    }
    // Hesseburg matrix after rotation
    m_Upper = Array<OneD, Array<OneD, NekDouble>>(m_LinSysMaxStorage);
    for (size_t nd = 0; nd < m_LinSysMaxStorage; nd++)
    {
        m_Upper[nd] = Array<OneD, NekDouble>(m_LinSysMaxStorage + 1, 0.0);
    }
    // Total search directions
    m_V_total = Array<OneD, Array<OneD, NekDouble>>(m_LinSysMaxStorage + 1);
}

References m_DifferenceFlag0, Nektar::LibUtilities::NekSysKey::m_DifferenceFlag0, m_DifferenceFlag1, Nektar::LibUtilities::NekSysKey::m_DifferenceFlag1, m_hes, m_KrylovMaxHessMatBand, Nektar::LibUtilities::NekSysKey::m_KrylovMaxHessMatBand, Nektar::LibUtilities::NekLinSysIter::m_LinSysMaxStorage, Nektar::LibUtilities::NekSys::m_maxiter, m_maxrestart, m_NekLinSysLeftPrecon, Nektar::LibUtilities::NekSysKey::m_NekLinSysLeftPrecon, m_NekLinSysRightPrecon, Nektar::LibUtilities::NekSysKey::m_NekLinSysRightPrecon, m_Upper, and m_V_total.

◆ ~NekLinSysIterGMRES()

Nektar::LibUtilities::NekLinSysIterGMRES::~NekLinSysIterGMRES ( )

overridedefault

Member Function Documentation

◆ create()

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

inlinestatic

Definition at line 52 of file NekLinSysIterGMRES.h.

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

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

◆ DoArnoldi()

void Nektar::LibUtilities::NekLinSysIterGMRES::DoArnoldi	(	const int	starttem,
		const int	endtem,
		const int	nGlobal,
		const int	nDir,
		Array< OneD, NekDouble > &	w,
		Array< OneD, NekDouble > &	Vsingle1,
		Array< OneD, NekDouble > &	Vsingle2,
		Array< OneD, NekDouble > &	hsingle
	)

private

Definition at line 422 of file NekLinSysIterGMRES.cpp.

{
    NekDouble alpha, beta, vExchange = 0.0;
    Array<OneD, NekDouble> tmp;
    int nNonDir = nGlobal - nDir;
    LibUtilities::Timer timer;
    timer.Start();
    m_operator.DoNekSysLhsEval(Vsingle1, w, m_DifferenceFlag1);
    timer.Stop();
    timer.AccumulateRegion("NekSysOperators::DoNekSysLhsEval", 10);
 
    if (m_NekLinSysLeftPrecon)
    {
        m_operator.DoNekSysPrecon(w + nDir, tmp = w + nDir);
    }
 
    Vmath::Smul(nNonDir, sqrt(m_prec_factor), w + nDir, 1, tmp = w + nDir, 1);
 
    // Modified Gram-Schmidt
    for (int i = starttem; i < endtem; ++i)
    {
        vExchange = Vmath::Dot2(nNonDir, &w[0] + nDir, &m_V_total[i][0] + nDir,
                                &m_map[0] + nDir);
        m_rowComm->AllReduce(vExchange, LibUtilities::ReduceSum);
 
        hsingle[i] = vExchange;
 
        beta = -1.0 * vExchange;
        Vmath::Svtvp(nNonDir, beta, &m_V_total[i][0] + nDir, 1, &w[0] + nDir, 1,
                     &w[0] + nDir, 1);
    }
    // end of Modified Gram-Schmidt
 
    // calculate the L2 norm and normalize
    vExchange =
        Vmath::Dot2(nNonDir, &w[0] + nDir, &w[0] + nDir, &m_map[0] + nDir);
    m_rowComm->AllReduce(vExchange, LibUtilities::ReduceSum);
 
    hsingle[endtem] = sqrt(vExchange);
 
    alpha = 1.0 / hsingle[endtem];
    Vmath::Smul(nNonDir, alpha, &w[0] + nDir, 1, &Vsingle2[0] + nDir, 1);
}

Referenced by DoGmresRestart().

◆ DoBackward()

void Nektar::LibUtilities::NekLinSysIterGMRES::DoBackward	(	const int	number,
		Array< OneD, Array< OneD, NekDouble > > &	A,
		const Array< OneD, const NekDouble > &	b,
		Array< OneD, NekDouble > &	y
	)

private

Definition at line 529 of file NekLinSysIterGMRES.cpp.

{
    // Number is the entry number
    // but C++'s order need to be one smaller
    int maxid = number - 1;
    NekDouble sum;
    y[maxid] = b[maxid] / A[maxid][maxid];
    for (int i = maxid - 1; i > -1; --i)
    {
        sum = b[i];
        for (int j = i + 1; j < number; ++j)
        {
            // i and j changes due to use Array<OneD,Array<OneD,NekDouble>>
            sum -= y[j] * A[j][i];
        }
        y[i] = sum / A[i][i];
    }
}

Referenced by DoGmresRestart().

◆ DoGivensRotation()

void Nektar::LibUtilities::NekLinSysIterGMRES::DoGivensRotation	(	const int	starttem,
		const int	endtem,
		const int	nGlobal,
		const int	nDir,
		Array< OneD, NekDouble > &	c,
		Array< OneD, NekDouble > &	s,
		Array< OneD, NekDouble > &	hsingle,
		Array< OneD, NekDouble > &	eta
	)

private

Definition at line 475 of file NekLinSysIterGMRES.cpp.

{
    NekDouble temp_dbl;
    NekDouble dd;
    NekDouble hh;
    int idtem = endtem - 1;
    // The starttem and endtem are beginning and ending order of Givens rotation
    // They usually equal to the beginning position and ending position of
    // Hessenburg matrix But sometimes starttem will change, like if it is
    // initial 0 and becomes nonzero because previous Givens rotation See Yu
    // Pan's User Guide
    for (int i = starttem; i < idtem; ++i)
    {
        temp_dbl       = c[i] * hsingle[i] - s[i] * hsingle[i + 1];
        hsingle[i + 1] = s[i] * hsingle[i] + c[i] * hsingle[i + 1];
        hsingle[i]     = temp_dbl;
    }
    dd = hsingle[idtem];
    hh = hsingle[endtem];
    if (hh == 0.0)
    {
        c[idtem] = 1.0;
        s[idtem] = 0.0;
    }
    else if (abs(hh) > abs(dd))
    {
        temp_dbl = -dd / hh;
        s[idtem] = 1.0 / sqrt(1.0 + temp_dbl * temp_dbl);
        c[idtem] = temp_dbl * s[idtem];
    }
    else
    {
        temp_dbl = -hh / dd;
        c[idtem] = 1.0 / sqrt(1.0 + temp_dbl * temp_dbl);
        s[idtem] = temp_dbl * c[idtem];
    }
 
    hsingle[idtem]  = c[idtem] * hsingle[idtem] - s[idtem] * hsingle[endtem];
    hsingle[endtem] = 0.0;
 
    temp_dbl    = c[idtem] * eta[idtem] - s[idtem] * eta[endtem];
    eta[endtem] = s[idtem] * eta[idtem] + c[idtem] * eta[endtem];
    eta[idtem]  = temp_dbl;
}

References tinysimd::abs(), and tinysimd::sqrt().

Referenced by DoGmresRestart().

◆ DoGMRES()

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

private

Actual iterative solve-GMRES.

Solve a global linear system using the Gmres 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.

Parameters

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

Definition at line 116 of file NekLinSysIterGMRES.cpp.

{
    m_prec_factor = NekConstants::kNekUnsetDouble;
 
    if (m_rhs_magnitude == NekConstants::kNekUnsetDouble)
    {
        NekVector<NekDouble> inGlob(nGlobal, pInput, eWrapper);
        Set_Rhs_Magnitude(inGlob);
    }
 
    // Get vector sizes
    NekDouble eps = 0.0;
    int nNonDir   = nGlobal - nDir;
 
    Array<OneD, NekDouble> tmp;
 
    // 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);
 
    m_totalIterations = 0;
    m_converged       = false;
 
    bool restarted = false;
    bool truncted  = false;
 
    if (m_KrylovMaxHessMatBand > 0)
    {
        truncted = true;
    }
 
    for (int nrestart = 0; nrestart < m_maxrestart; ++nrestart)
    {
        eps =
            DoGmresRestart(restarted, truncted, nGlobal, pInput, pOutput, nDir);
 
        if (m_converged)
        {
            break;
        }
        restarted = true;
    }
 
    if (m_verbose)
    {
        Array<OneD, NekDouble> r0(nGlobal, 0.0);
        m_operator.DoNekSysLhsEval(pOutput, r0, m_DifferenceFlag0);
        Vmath::Svtvp(nNonDir, -1.0, &r0[0] + nDir, 1, &pInput[0] + nDir, 1,
                     &r0[0] + nDir, 1);
        NekDouble vExchange = Vmath::Dot2(nNonDir, &r0[0] + nDir, &r0[0] + nDir,
                                          &m_map[0] + nDir);
        m_rowComm->AllReduce(vExchange, LibUtilities::ReduceSum);
        NekDouble eps1 = vExchange;
 
        if (m_root)
        {
            int nwidthcolm = 13;
 
            cout << std::scientific << std::setw(nwidthcolm)
                 << std::setprecision(nwidthcolm - 8)
                 << "       GMRES iterations made = " << m_totalIterations
                 << " using tolerance of " << m_tolerance
                 << " (error = " << sqrt(eps * m_prec_factor / m_rhs_magnitude)
                 << ")";
 
            cout << " WITH (GMRES eps = " << eps << " REAL eps= " << eps1
                 << ")";
 
            if (m_converged)
            {
                cout << " CONVERGED" << endl;
            }
            else
            {
                cout << " WARNING: Exceeded maxIt" << endl;
            }
        }
    }
 
    if (m_FlagWarnings)
    {
        WARNINGL1(m_converged, "GMRES did not converge.");
    }
    return m_totalIterations;
}

Referenced by v_SolveSystem().

◆ DoGmresRestart()

NekDouble Nektar::LibUtilities::NekLinSysIterGMRES::DoGmresRestart	(	const bool	restarted,
		const bool	truncted,
		const int	nGlobal,
		const Array< OneD, const NekDouble > &	pInput,
		Array< OneD, NekDouble > &	pOutput,
		const int	nDir
	)

private

Actual iterative gmres solver for one restart.

Definition at line 205 of file NekLinSysIterGMRES.cpp.

{
    int nNonDir = nGlobal - nDir;
 
    // Allocate array storage of coefficients
    // Residual
    Array<OneD, NekDouble> eta(m_LinSysMaxStorage + 1, 0.0);
    // Givens rotation c
    Array<OneD, NekDouble> cs(m_LinSysMaxStorage, 0.0);
    // Givens rotation s
    Array<OneD, NekDouble> sn(m_LinSysMaxStorage, 0.0);
    // Total coefficients, just for check
    Array<OneD, NekDouble> y_total(m_LinSysMaxStorage, 0.0);
    // Search direction order
    Array<OneD, int> id(m_LinSysMaxStorage, 0);
    Array<OneD, int> id_start(m_LinSysMaxStorage, 0);
    Array<OneD, int> id_end(m_LinSysMaxStorage, 0);
    // Temporary variables
    int idtem;
    int starttem;
    int endtem;
 
    NekDouble eps;
    NekDouble beta, alpha;
    NekDouble vExchange = 0;
    // Temporary Array
    Array<OneD, NekDouble> w(nGlobal, 0.0);
    Array<OneD, NekDouble> wk(nGlobal, 0.0);
    Array<OneD, NekDouble> r0(nGlobal, 0.0);
    Array<OneD, NekDouble> tmp;
    Array<OneD, NekDouble> Vsingle1;
    Array<OneD, NekDouble> Vsingle2;
    Array<OneD, NekDouble> hsingle1;
    Array<OneD, NekDouble> hsingle2;
 
    if (restarted)
    {
        // This is A*x
        m_operator.DoNekSysLhsEval(pOutput, r0, m_DifferenceFlag0);
 
        // The first search direction
        beta = -1.0;
 
        // This is r0 = b-A*x
        Vmath::Svtvp(nNonDir, beta, &r0[0] + nDir, 1, &pInput[0] + nDir, 1,
                     &r0[0] + nDir, 1);
    }
    else
    {
        // This is r0 = b, x = x0 assumed to be zero
        Vmath::Vcopy(nNonDir, &pInput[0] + nDir, 1, &r0[0] + nDir, 1);
    }
 
    if (m_NekLinSysLeftPrecon)
    {
        m_operator.DoNekSysPrecon(r0 + nDir, tmp = r0 + nDir);
    }
 
    // Norm of (r0)
    // The m_map tells how to connect
    vExchange =
        Vmath::Dot2(nNonDir, &r0[0] + nDir, &r0[0] + nDir, &m_map[0] + nDir);
    m_rowComm->AllReduce(vExchange, LibUtilities::ReduceSum);
    eps = vExchange;
 
    // Detect zero input array
    // Causes Arnoldi to breakdown, hence stop here
    if (eps < m_tolerance * m_tolerance * m_rhs_magnitude)
    {
        m_converged = true;
        if (m_prec_factor == NekConstants::kNekUnsetDouble)
        {
            m_prec_factor = 1.0;
        }
        return eps;
    }
 
    if (!restarted)
    {
        if (m_NekLinSysLeftPrecon &&
            m_prec_factor == NekConstants::kNekUnsetDouble)
        {
            // Evaluate initial residual error for exit check
            vExchange = Vmath::Dot2(nNonDir, &pInput[0] + nDir,
                                    &pInput[0] + nDir, &m_map[0] + nDir);
            m_rowComm->AllReduce(vExchange, LibUtilities::ReduceSum);
            m_prec_factor = vExchange / eps;
        }
        else
        {
            m_prec_factor = 1.0;
        }
    }
 
    Vmath::Smul(nNonDir, sqrt(m_prec_factor), r0 + nDir, 1, tmp = r0 + nDir, 1);
    eps    = eps * m_prec_factor;
    eta[0] = sqrt(eps);
 
    // Give an order for the entries in Hessenburg matrix
    for (int nd = 0; nd < m_LinSysMaxStorage; ++nd)
    {
        id[nd]     = nd;
        id_end[nd] = nd + 1;
        starttem   = id_end[nd] - m_KrylovMaxHessMatBand;
        if (truncted && (starttem) > 0)
        {
            id_start[nd] = starttem;
        }
        else
        {
            id_start[nd] = 0;
        }
    }
 
    // Normlized by r0 norm V(:,1)=r0/norm(r0)
    alpha = 1.0 / eta[0];
 
    // Scalar multiplication
    if (m_V_total[0].size() == 0)
    {
        m_V_total[0] = Array<OneD, NekDouble>(nGlobal, 0.0);
    }
    Vmath::Smul(nNonDir, alpha, &r0[0] + nDir, 1, &m_V_total[0][0] + nDir, 1);
 
    // Restarted Gmres(m) process
    if (m_NekLinSysRightPrecon)
    {
        Vsingle1 = Array<OneD, NekDouble>(nGlobal, 0.0);
    }
 
    int nswp = 0;
    for (int nd = 0; nd < m_LinSysMaxStorage; ++nd)
    {
        if (m_V_total[nd + 1].size() == 0)
        {
            m_V_total[nd + 1] = Array<OneD, NekDouble>(nGlobal, 0.0);
        }
        Vmath::Zero(nGlobal, m_V_total[nd + 1], 1);
        Vmath::Zero(m_LinSysMaxStorage + 1, m_hes[nd], 1);
        Vsingle2 = m_V_total[nd + 1];
        hsingle1 = m_hes[nd];
 
        if (m_NekLinSysRightPrecon)
        {
            m_operator.DoNekSysPrecon(m_V_total[nd] + nDir,
                                      tmp = Vsingle1 + nDir);
        }
        else
        {
            Vsingle1 = m_V_total[nd];
        }
 
        // w here is no need to add nDir due to temporary Array
        idtem    = id[nd];
        starttem = id_start[idtem];
        endtem   = id_end[idtem];
 
        DoArnoldi(starttem, endtem, nGlobal, nDir, w, Vsingle1, Vsingle2,
                  hsingle1);
 
        if (starttem > 0)
        {
            starttem = starttem - 1;
        }
 
        hsingle2 = m_Upper[nd];
        Vmath::Vcopy(m_LinSysMaxStorage + 1, &hsingle1[0], 1, &hsingle2[0], 1);
        DoGivensRotation(starttem, endtem, nGlobal, nDir, cs, sn, hsingle2,
                         eta);
 
        eps = eta[nd + 1] * eta[nd + 1];
 
        // This Gmres merge truncted Gmres to accelerate.
        // If truncted, cannot jump out because
        // the last term of eta is not residual
        if ((!truncted) || (nd < m_KrylovMaxHessMatBand))
        {
            if ((eps < m_tolerance * m_tolerance * m_rhs_magnitude) && nd > 0)
            {
                m_converged = true;
            }
        }
        nswp++;
        m_totalIterations++;
 
        if (m_converged)
        {
            break;
        }
    }
 
    DoBackward(nswp, m_Upper, eta, y_total);
 
    // Calculate output V_total * y_total.
    Array<OneD, NekDouble> solution(nNonDir, 0.0);
    for (int i = 0; i < nswp; ++i)
    {
        Vmath::Svtvp(nNonDir, y_total[i], &m_V_total[i][0] + nDir, 1,
                     solution.get(), 1, solution.get(), 1);
    }
 
    if (m_NekLinSysRightPrecon)
    {
        m_operator.DoNekSysPrecon(solution, solution);
    }
 
    // Update output.
    Vmath::Vadd(nNonDir, solution.get(), 1, &pOutput[0] + nDir, 1,
                &pOutput[0] + nDir, 1);
 
    return eps;
}

Referenced by DoGMRES().

◆ GetMaxLinIte()

int Nektar::LibUtilities::NekLinSysIterGMRES::GetMaxLinIte ( )

inline

Definition at line 72 of file NekLinSysIterGMRES.h.

    {
        return (m_maxrestart * m_LinSysMaxStorage);
    }

References Nektar::LibUtilities::NekLinSysIter::m_LinSysMaxStorage, and m_maxrestart.

◆ v_InitObject()

void Nektar::LibUtilities::NekLinSysIterGMRES::v_InitObject ( )

overrideprotectedvirtual

Reimplemented from Nektar::LibUtilities::NekLinSysIter.

Definition at line 86 of file NekLinSysIterGMRES.cpp.

{
    NekLinSysIter::v_InitObject();
}

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

◆ v_SolveSystem()

int Nektar::LibUtilities::NekLinSysIterGMRES::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 94 of file NekLinSysIterGMRES.cpp.

{
    m_tolerance     = max(tol, 1.0E-15);
    int niterations = DoGMRES(nGlobal, pInput, pOutput, nDir);
 
    return niterations;
}

References DoGMRES(), and Nektar::LibUtilities::NekSys::m_tolerance.

Friends And Related Function Documentation

◆ MemoryManager< NekLinSysIterGMRES >

friend class MemoryManager< NekLinSysIterGMRES >

friend

Support creation through MemoryManager.

Definition at line 1 of file NekLinSysIterGMRES.h.

Member Data Documentation

◆ className

string Nektar::LibUtilities::NekLinSysIterGMRES::className

static

Initial value:

=
    LibUtilities::GetNekLinSysIterFactory().RegisterCreatorFunction(
        "GMRES", NekLinSysIterGMRES::create, "NekLinSysIterGMRES solver.")

Definition at line 64 of file NekLinSysIterGMRES.h.

◆ def

std::string Nektar::LibUtilities::NekLinSysIterGMRES::def

staticprivate

Definition at line 130 of file NekLinSysIterGMRES.h.

◆ lookupIds

std::string Nektar::LibUtilities::NekLinSysIterGMRES::lookupIds[]

staticprivate

Definition at line 129 of file NekLinSysIterGMRES.h.

◆ m_DifferenceFlag0

bool Nektar::LibUtilities::NekLinSysIterGMRES::m_DifferenceFlag0 = false

protected

Definition at line 88 of file NekLinSysIterGMRES.h.

Referenced by DoGMRES(), DoGmresRestart(), and NekLinSysIterGMRES().

◆ m_DifferenceFlag1

bool Nektar::LibUtilities::NekLinSysIterGMRES::m_DifferenceFlag1 = false

protected

Definition at line 89 of file NekLinSysIterGMRES.h.

Referenced by DoArnoldi(), and NekLinSysIterGMRES().

◆ m_hes

Array<OneD, Array<OneD, NekDouble> > Nektar::LibUtilities::NekLinSysIterGMRES::m_hes

private

Definition at line 133 of file NekLinSysIterGMRES.h.

Referenced by DoGmresRestart(), and NekLinSysIterGMRES().

◆ m_KrylovMaxHessMatBand

int Nektar::LibUtilities::NekLinSysIterGMRES::m_KrylovMaxHessMatBand

protected

Definition at line 83 of file NekLinSysIterGMRES.h.

Referenced by DoGMRES(), DoGmresRestart(), and NekLinSysIterGMRES().

◆ m_maxrestart

int Nektar::LibUtilities::NekLinSysIterGMRES::m_maxrestart

protected

Definition at line 79 of file NekLinSysIterGMRES.h.

Referenced by DoGMRES(), GetMaxLinIte(), and NekLinSysIterGMRES().

◆ m_NekLinSysLeftPrecon

bool Nektar::LibUtilities::NekLinSysIterGMRES::m_NekLinSysLeftPrecon = false

protected

Definition at line 85 of file NekLinSysIterGMRES.h.

Referenced by DoArnoldi(), DoGmresRestart(), and NekLinSysIterGMRES().

◆ m_NekLinSysRightPrecon

bool Nektar::LibUtilities::NekLinSysIterGMRES::m_NekLinSysRightPrecon = true

protected

Definition at line 86 of file NekLinSysIterGMRES.h.

Referenced by DoGmresRestart(), and NekLinSysIterGMRES().

◆ m_Upper

Array<OneD, Array<OneD, NekDouble> > Nektar::LibUtilities::NekLinSysIterGMRES::m_Upper

private

Definition at line 135 of file NekLinSysIterGMRES.h.

Referenced by DoGmresRestart(), and NekLinSysIterGMRES().

◆ m_V_total

Array<OneD, Array<OneD, NekDouble> > Nektar::LibUtilities::NekLinSysIterGMRES::m_V_total

private

Definition at line 137 of file NekLinSysIterGMRES.h.

Referenced by DoArnoldi(), DoGmresRestart(), and NekLinSysIterGMRES().

Public Member Functions

Static Public Member Functions

Static Public Attributes

Protected Member Functions

Protected Attributes

Private Member Functions

Private Attributes

Static Private Attributes

Friends

Detailed Description

Constructor & Destructor Documentation

◆ NekLinSysIterGMRES()

◆ ~NekLinSysIterGMRES()

Member Function Documentation

◆ create()

◆ DoArnoldi()

◆ DoBackward()

◆ DoGivensRotation()

◆ DoGMRES()

◆ DoGmresRestart()

◆ GetMaxLinIte()

◆ v_InitObject()

◆ v_SolveSystem()

Friends And Related Function Documentation

◆ MemoryManager< NekLinSysIterGMRES >

Member Data Documentation

◆ className

◆ def

◆ lookupIds

◆ m_DifferenceFlag0

◆ m_DifferenceFlag1

◆ m_hes

◆ m_KrylovMaxHessMatBand

◆ m_maxrestart

◆ m_NekLinSysLeftPrecon

◆ m_NekLinSysRightPrecon

◆ m_Upper

◆ m_V_total