dg::ChebyshevIteration< ContainerType > Class Template Reference

Preconditioned Chebyshev iteration for solving \( PAx=Pb\). More...

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

Public Member Functions
	ChebyshevIteration ()=default
	Allocate nothing, Call construct method before usage. More...
	ChebyshevIteration (const ContainerType &copyable)
	Allocate memory for the pcg method. More...
const ContainerType &	copyable () const
	Return an object of same size as the object used for construction. More...
void	construct (const ContainerType &copyable)
	Allocate memory for the pcg method. More...
template<class MatrixType , class ContainerType0 , class ContainerType1 >
void	solve (MatrixType &&A, ContainerType0 &x, const ContainerType1 &b, value_type min_ev, value_type max_ev, unsigned num_iter, bool x_is_zero=false)
	Solve the system \( Ax = b\) using num_iter Chebyshev iteration. More...
template<class MatrixType0 , class MatrixType1 , class ContainerType0 , class ContainerType1 >
void	solve (MatrixType0 &&A, ContainerType0 &x, const ContainerType1 &b, MatrixType1 &&P, value_type min_ev, value_type max_ev, unsigned num_iter, bool x_is_zero=false)
	Solve the system \( PAx = Pb\) using num_iter Preconditioned Chebyshev iteration. More...

Detailed Description

template<class ContainerType>
class dg::ChebyshevIteration< ContainerType >

Preconditioned Chebyshev iteration for solving \( PAx=Pb\).

Chebyshev iteration is not well-suited for solving matrix equations on its own. Rather, it is suited as a smoother for a multigrid algorithm and also as a Preconditioner for the Conjugate Gradient method. It does not contain scalar products, which makes it appaeling for both small and highly parallelized systems.

Given the minimum and maximum Eigenvalue of the matrix A we define

\[ \theta = (\lambda_\min+\lambda_\max)/2 \quad \delta = (\lambda_\max - \lambda_\min)/2 \\ \rho_0 := \frac{\delta}{\theta},\ x_0 := x, \ x_{1} = x_0+\frac{1}{\theta} P(b-Ax_0) \\ \rho_{k}:=\left(\frac{2\theta}{\delta}-\rho_{k-1}\right)^{-1} \\ x_{k+1} := x_k + \rho_k\left( \rho_{k-1}(x_k - x_{k-1}) + \frac{2}{\delta} P( b - Ax_k) \right) \]

The preconditioned version is obtained by applying the regular version to \( \bar A\bar x = \bar b\) with \( \bar A := {E^{-1}}^\mathrm{T} A E^{-1} \), \( \bar x := Ex\) and \( \bar b := {E^{-1}}^\mathrm{T}\), where \( P = {E^{-1}}^\mathrm{T} E^{-1}\) is the preconditioner. The bounds on the spectrum then need to be on the \(PA\) matrix.

Note

The maximum Eigenvalue of \( A\) and \( P A\) can be estimated using the EVE class.

See also

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

Note

If the initial vector is zero Chebyshev iteration will produce the Chebyshev polynomial \( C_k( A) b\) applied to the right hand side and the preconditioned version produces \( C_{k-1}(PA)Pb = E^{-1} C_{k-1}( {E^{-1}}^\mathrm{T} A E^{-1}){E^{-1}}^\mathrm{T}\)

Attention

Chebyshev iteration may diverge if the elliptical bound of the Eigenvalues is not accurate (in particular if \(\lambda_\max\) is underestimated) or if an ellipsis is not a good fit for the spectrum of the matrix

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::ChebyshevIteration< ContainerType >::container_type = ContainerType

value_type

template<class ContainerType >

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

value type of the ContainerType class

Constructor & Destructor Documentation

ChebyshevIteration() [1/2]

template<class ContainerType >

dg::ChebyshevIteration< ContainerType >::ChebyshevIteration ( )

default

Allocate nothing, Call construct method before usage.

ChebyshevIteration() [2/2]

template<class ContainerType >

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

inline

Allocate memory for the pcg method.

Parameters

copyable

A ContainerType must be copy-constructible from this

Member Function Documentation

construct()

template<class ContainerType >

void dg::ChebyshevIteration< ContainerType >::construct ( const ContainerType &  copyable )

inline

Allocate memory for the pcg method.

Parameters

copyable

A ContainerType must be copy-constructible from this

copyable()

template<class ContainerType >

const ContainerType & dg::ChebyshevIteration< 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

solve() [1/2]

template<class ContainerType >

template<class MatrixType , class ContainerType0 , class ContainerType1 >

void dg::ChebyshevIteration< ContainerType >::solve ( MatrixType &&  A, ContainerType0 &  x, const ContainerType1 &  b, value_type  min_ev, value_type  max_ev, unsigned  num_iter, bool  x_is_zero = false  )

inline

Solve the system \( Ax = b\) using num_iter Chebyshev iteration.

The iteration stops when the maximum number of iterations is reached

Parameters

A	A symmetric, positive definit 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.
min_ev	an estimate of the minimum Eigenvalue
max_ev	an estimate of the maximum Eigenvalue of \( A\) (must be larger than min_ev) Use EVE to get this value
num_iter	the number of iterations k (equals the number of times A is applied) If 0 the function returns immediately
x_is_zero	If true, the first matrix-vector multiplication is avoided by assuming x is zero. (This works even if x is not actually 0) This is in particular in the case when Chebyshev Iteration is used as a Preconditioner

Note

In the x_is_zero mode k iterations will produce the k-1 Chebyshev polynomial applied to the right hand side \( x = C_{k-1}(A)b\)

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

solve() [2/2]

template<class ContainerType >

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

void dg::ChebyshevIteration< ContainerType >::solve ( MatrixType0 &&  A, ContainerType0 &  x, const ContainerType1 &  b, MatrixType1 &&  P, value_type  min_ev, value_type  max_ev, unsigned  num_iter, bool  x_is_zero = false  )

inline

Solve the system \( PAx = Pb\) using num_iter Preconditioned Chebyshev iteration.

The iteration stops when the maximum number of iterations is reached

Parameters

A	A symmetric, positive definit 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 ( \( M^{-1}\) in the above notation
min_ev	an estimate of the minimum Eigenvalue
max_ev	an estimate of the maximum Eigenvalue of \( P A\) (must be larger than min_ev) Use EVE to get this value
num_iter	the number of iterations k (equals the number of times A is applied) If 0 the function returns immediately
x_is_zero	If true, the first matrix-vector multiplication is avoided by assuming x is zero. (This works even if x is not actually 0) This is in particular in the case when Chebyshev Iteration is used as a Preconditioner

Note

In the x_is_zero mode k iterations will produce the k-1 Chebyshev polynomial applied to the right hand side \( x = C_{k-1}(PA)Pb = E^{-1} C_{k-1}( {E^{-1}}^\mathrm{T} A E^{-1}){E^{-1}}^\mathrm{T}\)

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:

• chebyshev.h