|
| | LGMRES ()=default |
| | Allocate nothing, Call construct method before usage.
|
| |
| | LGMRES (const ContainerType ©able, unsigned max_inner, unsigned max_outer, unsigned max_restarts) |
| | Allocate memory for the preconditioned LGMRES method.
|
| |
| template<class ... Params> |
| void | construct (Params &&...ps) |
| | Perfect forward parameters to one of the constructors.
|
| |
| void | set_max (unsigned new_Restarts) |
| | Set the number of restarts.
|
| |
| unsigned | get_max () const |
| | Get the current maximum number of restarts.
|
| |
| void | set_throw_on_fail (bool throw_on_fail) |
| | Set or unset a throw on failure-to-converge.
|
| |
| const ContainerType & | copyable () const |
| | Return an object of same size as the object used for construction.
|
| |
| 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.
|
| |
| bool | converged () const |
| | If last call to solve converged or not.
|
| |
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
unsigned num_iter = lgmres.solve( A, x, b, A.precond(), A.weights(), 1e-6);
INFO( "After "<<num_iter<<" LGMRES iterations we have");
INFO( "L2 Norm of Error is " << res);
CHECK( res < 1e-6);
- 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
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
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. |
- 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
1.12.0