dg::PCG< ContainerType > Class Template Reference

Preconditioned conjugate gradient method to solve \( Ax=b\). More...

Public Types
using	container_type = ContainerType
using	value_type = get_value_type< ContainerType >

Public Member Functions
	PCG ()=default
	Allocate nothing, Call construct method before usage. More...
	PCG (const ContainerType &copyable, unsigned max_iterations)
	Allocate memory for the pcg method. More...
void	set_max (unsigned new_max)
	Set the maximum number of iterations. More...
unsigned	get_max () const
	Get the current maximum number of iterations. More...
const ContainerType &	copyable () const
	Return an object of same size as the object used for construction. More...
void	set_verbose (bool verbose)
	Set or unset debugging output during iterations. More...
void	set_throw_on_fail (bool throw_on_fail)
	Set or unset a throw on failure-to-converge. More...
template<class ... Params>
void	construct (Params &&...ps)
	Perfect forward parameters to one of the constructors. 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, int test_frequency=1)
	Solve \( Ax = b\) using a preconditioned conjugate gradient method. More...

Detailed Description

template<class ContainerType>
class dg::PCG< ContainerType >

Preconditioned conjugate gradient method to solve \( Ax=b\).

where \( A\) is positive definite and self-adjoint in the weighted scalar product (defined by the diagonal weights matrix \(W\))

\[ A^\dagger := \frac{1}{W} A^T W = A\]

. Note that if \( A\) is self-adjoint then both \( (WA)^T = WA \) and \( \left(A \frac{1}{W}\right)^T = A\frac{1}{W}\) are symmetric. The positive definite, self-adjoint preconditioner \( P \approx A^{-1}\) that approximates the inverse of \( A\) and is fast to apply, is used to solve the left preconditioned system

\[ PAx=Pb\]

Note

Our implementation uses a stopping criterion based on the residual at iteration i \( || r_i ||_W = ||Ax_i-b||_W < \epsilon( ||b||_W + C) \). However, the real error is bound by

\[ \frac{ ||e_i||_W}{||x||_W} \leq \kappa(PA) \frac{||r_i||_W}{||b||_W} \]

Thus, if the condition number \( \kappa\) is large the real error may still be large even if the residual error is small see Ashby et al., The Role of the Inner Product in Stopping Criteria for Conjugate Gradient Iterations (2001)

See also

This implements the PCG algorithm (applied to \((WA)\) as given in https://en.wikipedia.org/wiki/Conjugate_gradient_method or the book Iterative Methods for Sparse Linear Systems" 2nd edition by Yousef Saad

Note

Conjugate gradients might become unstable for positive semidefinite matrices arising e.g. in the discretization of the periodic laplacian

Attention

beware the sign: a negative definite matrix does not work in Conjugate gradient

    // create volume and inverse volume on previously defined grid
    const dg::HVec w2d = dg::create::weights( grid);
 
    // Create unnormalized Laplacian
    dg::Elliptic<dg::CartesianGrid2d, dg::HMatrix, dg::HVec> A( grid);
 
    // allocate memory in conjugate gradient
    dg::PCG<dg::HVec > pcg( copyable_vector, max_iter);
 
    // Evaluate right hand side and solution on the grid
    dg::HVec b = dg::evaluate ( laplace_fct, grid);
    const dg::HVec solution = dg::evaluate ( fct, grid);
 
    // use inverse volume as preconditioner in solution method
    const double eps = 1e-6;
    unsigned num_iter = pcg.solve( A, x, b, 1., w2d, eps);

Template Parameters

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

Member Typedef Documentation

container_type

template<class ContainerType >

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

value_type

template<class ContainerType >

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

value type of the ContainerType class

Constructor & Destructor Documentation

PCG() [1/2]

template<class ContainerType >

dg::PCG< ContainerType >::PCG ( )

default

Allocate nothing, Call construct method before usage.

PCG() [2/2]

template<class ContainerType >

dg::PCG< ContainerType >::PCG ( const ContainerType &  copyable, unsigned  max_iterations  )

inline

Allocate memory for the pcg method.

Parameters

copyable	A ContainerType must be copy-constructible from this
max_iterations	Maximum number of iterations to be used

Member Function Documentation

construct()

template<class ContainerType >

template<class ... Params>

void dg::PCG< 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

copyable()

template<class ContainerType >

const ContainerType & dg::PCG< 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::PCG< ContainerType >::get_max ( ) const

inline

Get the current maximum number of iterations.

Returns

the current maximum

set_max()

template<class ContainerType >

void dg::PCG< ContainerType >::set_max ( unsigned  new_max )

inline

Set the maximum number of iterations.

Parameters

new_max

New maximum number

set_throw_on_fail()

template<class ContainerType >

void dg::PCG< 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

set_verbose()

template<class ContainerType >

void dg::PCG< ContainerType >::set_verbose ( bool  verbose )

inline

Set or unset debugging output during iterations.

Parameters

verbose

If true, additional output will be written to std::cout during solution

solve()

template<class ContainerType >

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

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

Solve \( Ax = b\) using a preconditioned conjugate gradient 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 self-adjoint positive definit matrix with respect to the weights W
x	Contains an initial value on input and the solution on output.
b	The right hand side vector.
P	The preconditioner to be used (an approximation to the inverse of A that is fast to compute)
W	Weights that define the scalar product in which A and P are self-adjoint and in which the error norm is computed.
eps	The relative error to be respected
nrmb_correction	the absolute error C in units of eps to be respected
test_frequency	if set to 1 then the norm of the error is computed in every iteration to test if the loop can be terminated. Sometimes, especially for small sizes the dot product is expensive to compute, then it is beneficial to set this parameter to e.g. 10, which means that the errror condition is only evaluated every 10th iteration.

Returns

Number of iterations used to achieve desired precision

Note

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

Required memops per iteration (P is assumed vector):

• 15 reads + 4 writes
• plus the number of memops for A;

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:

• pcg.h