Functor class for the right preconditioned LGMRES method to solve \( Ax=b\). More...

Public Types
using	container_type = ContainerType

using	value_type = dg::get_value_type< ContainerType >

Public Member Functions
	LGMRES ()=default
	Allocate nothing, Call `construct` method before usage. More...

	LGMRES (const ContainerType &copyable, unsigned max_inner, unsigned max_outer, unsigned max_restarts)
	Allocate memory for the preconditioned LGMRES method. More...

template<class ... Params>
void	construct (Params &&...ps)
	Perfect forward parameters to one of the constructors. More...

void	set_max (unsigned new_Restarts)
	Set the number of restarts. More...

unsigned	get_max () const
	Get the current maximum number of restarts. More...

void	set_throw_on_fail (bool throw_on_fail)
	Set or unset a throw on failure-to-converge. More...

const ContainerType &	copyable () const
	Return an object of same size as the object used for construction. More...

template<class MatrixType0 , class ContainerType0 , class ContainerType1 , class MatrixType1 , class ContainerType2 >
unsigned	solve (MatrixType0 &&A, ContainerType0 &x, const ContainerType1 &b, MatrixType1 &&P, const ContainerType2 &W, value_type eps=1e-12, value_type nrmb_correction=1)
	Solve \( Ax = b\) using a right preconditioned LGMRES method. More...

bool	converged () const
	If last call to solve converged or not. More...

Detailed Description

template<class ContainerType>
class dg::LGMRES< ContainerType >

Functor class for the right preconditioned LGMRES method to solve \( Ax=b\).

Note: GMRES is a method for solving non-symmetrical linear systems. LGMRES is a modification of restarted GMRES that aims to improve convergence. This implementation is adapted from: https://github.com/KellyBlack/GMRES/blob/master/GMRES.h with LGMRES elements from https://github.com/haranjackson/NewtonKrylov/blob/master/lgmres.cpp A paper can be found at https://www.cs.colorado.edu/~jessup/SUBPAGES/PS/lgmres.pdf

See also: For more information see the book Iteratvie Methods for Sparse Linear Systems" 2nd edition by Yousef Saad

Basically the Krylov subspace is augmented by k approximations to the error. The storage requirement is m+3k vectors

Note: We use right preconditioning because this makes the residual norm automatically available in each iteration; the orthogonalization is done with respect to a user-provided inner product W.; the first cycle of LGMRES(m,k) is equivalent to the first cycle of GMRES(m+k); Only m matrix-vector multiplications need to be computed in LGMRES(m,k) per restart cycle irrespective of the value of k

Attention: There are a lot of calls to dot in LGMRES such that GPUs may struggle for small vector sizes; LGMRES can stagnate if the matrix A is not positive definite. The use of a preconditioner is paramount in such a situation

Member Typedef Documentation

◆ container_type

template<class ContainerType >

using dg::LGMRES< ContainerType >::container_type = ContainerType

◆ value_type

template<class ContainerType >

using dg::LGMRES< ContainerType >::value_type = dg::get_value_type<ContainerType>

value type of the ContainerType class

Constructor & Destructor Documentation

◆ LGMRES() [1/2]

template<class ContainerType >

dg::LGMRES< ContainerType >::LGMRES ( )

default

Allocate nothing, Call construct method before usage.

◆ LGMRES() [2/2]

template<class ContainerType >

dg::LGMRES< ContainerType >::LGMRES	(	const ContainerType &	copyable,
		unsigned	max_inner,
		unsigned	max_outer,
		unsigned	max_restarts
	)

inline

Allocate memory for the preconditioned LGMRES method.

Parameters

copyable	A ContainerType must be copy-constructible from this
max_inner	Maximum number inner gmres iterations per restart. Usually 20-30 seems to be a decent number. Per iteration a matrix-vector product and a preconditioner-vector product needs to be computed.
max_outer	Maximum number of solutions (actually approximations to the error) saved for restart. Usually 1...3 is a good number. The Krylov Dimension is thus augmented to `max_inner+max_outer`. No new matrix-vector products need to be computed for the additional solutions. `max_outer=0` corresponds to standard GMRES.
max_restarts	Maximum number of restarts. The total maximum number of iterations/ matrix-vector products is thus `max_restarts*max_inner`

Member Function Documentation

◆ construct()

template<class ContainerType >

template<class ... Params>

void dg::LGMRES< ContainerType >::construct ( Params &&... ps )

inline

Perfect forward parameters to one of the constructors.

Template Parameters

Params deduced by the compiler

Parameters

ps	parameters forwarded to constructors

◆ converged()

template<class ContainerType >

bool dg::LGMRES< ContainerType >::converged ( ) const

inline

If last call to solve converged or not.

Returns: true if convergence was reached, false else

◆ copyable()

template<class ContainerType >

const ContainerType & dg::LGMRES< ContainerType >::copyable ( ) const

inline

Return an object of same size as the object used for construction.

Returns: A copyable object; what it contains is undefined, its size is important

◆ get_max()

template<class ContainerType >

unsigned dg::LGMRES< ContainerType >::get_max ( ) const

inline

Get the current maximum number of restarts.

Returns: the current maximum of restarts

◆ set_max()

template<class ContainerType >

void dg::LGMRES< ContainerType >::set_max ( unsigned new_Restarts )

inline

Set the number of restarts.

Parameters

new_Restarts New maximum number of restarts

◆ set_throw_on_fail()

template<class ContainerType >

void dg::LGMRES< ContainerType >::set_throw_on_fail ( bool throw_on_fail )

inline

Set or unset a throw on failure-to-converge.

Parameters

throw_on_fail If true, the solve method will thow a dg::Fail if it is unable to converge

Note: the default value is true

◆ solve()

template<class ContainerType >

template<class MatrixType0 , class ContainerType0 , class ContainerType1 , class MatrixType1 , class ContainerType2 >

unsigned dg::LGMRES< ContainerType >::solve	(	MatrixType0 &&	A,
		ContainerType0 &	x,
		const ContainerType1 &	b,
		MatrixType1 &&	P,
		const ContainerType2 &	W,
		value_type	eps = `1e-12`,
		value_type	nrmb_correction = `1`
	)

Solve \( Ax = b\) using a right preconditioned LGMRES method.

The iteration stops if \( ||Ax-b||_W < \epsilon( ||b||_W + C) \) where \(C\) is the absolute error in units of \( \epsilon\) and \( W \) defines a square norm

Parameters

A	A matrix
x	Contains an initial value on input and the solution on output.
b	The right hand side vector. x and b may be the same vector.
P	The preconditioner to be used
W	A diagonal matrix (a vector) that is used to define the scalar product in which the orthogonalization in LGMRES is computed and that defines the norm in which the stopping criterion is computed
eps	The relative error to be respected
nrmb_correction	the absolute error `C` in units of `eps` to be respected

Returns: Number of times the matrix A and the preconditioner P were multiplied to achieve the desired precision

Note: The method will throw dg::Fail if the desired accuracy is not reached within max_restarts You can unset this behaviour with the set_throw_on_fail member

Template Parameters

MatrixType

Any class for which a specialization of TensorTraits exists and defines a tensor_category derived from AnyMatrixTag. Furthermore, any functor/lambda type with signature void operator()( const ContainerType0&, ContainerType1&) . For example

scalar or Container: Scalars and containers act as diagonal matrices.
std::function<void( const ContainerType0&, ContainerType1&)>
dg::HMatrix and dg::IHMatrix with dg::HVec or std::vector<dg::HVec>
dg::DMatrix and dg::IDMatrix with dg::DVec or std::vector<dg::DVec>
dg::MHMatrix with dg::MHVec or std::vector<dg::MHVec>
dg::MDMatrix with dg::MDVec or std::vector<dg::MDVec>

In case of SelfMadeMatrixTag only those blas2 functions that have a corresponding member function in the Matrix class (e.g. symv( const ContainerType0&, ContainerType1&); ) can be called. If a Container has the RecursiveVectorTag, then the matrix is applied to each of the elements unless the type has the SelfMadeMatrixTag or is a Functor type.

ContainerType

Any class for which a specialization of TensorTraits exists and which fulfills the requirements of the there defined data and execution policies derived from AnyVectorTag and AnyPolicyTag. Among others

dg::HVec (serial), dg::DVec (cuda / omp), dg::MHVec (mpi + serial) or dg::MDVec (mpi + cuda / omp)
std::vector<dg::DVec> (vector of shared device vectors), std::array<double, 4> (array of 4 doubles) or std::map < std::string, dg::DVec> ( a map of named vectors)
double (scalar) and other primitive types ...

If there are several ContainerTypes in the argument list, then TensorTraits must exist for all of them

See also: See The dg dispatch system for a detailed explanation of our type dispatch system

The documentation for this class was generated from the following file:

lgmres.h

Public Types

Public Member Functions

Detailed Description

Member Typedef Documentation

◆ container_type

◆ value_type

Constructor & Destructor Documentation

◆ LGMRES() [1/2]

◆ LGMRES() [2/2]

Member Function Documentation

◆ construct()

◆ converged()

◆ copyable()

◆ get_max()

◆ set_max()

◆ set_throw_on_fail()

◆ solve()