Iterative Solvers

An iterative solver is a Krylov subspace method that may be combined with a preconditioner in order to solve the linear system.

All iterative solvers in AMGCL have two template parameters, Backend and InnerProduct. The Backend template parameter specifies the backend to target, and the InnerProduct parameter is used to select the implementation of the inner product to use with the solver. The correct implementation should be automatically selected by the library depending on whether the solver is used in a shared or distributed memory setting.

All solvers provide similar interface described below:

constructor(size_t n, const params &prm = params(), const backend_params &bprm = backend_params())

The solver constructor. Takes the size of the system to solve, the solver parameters and the backend parameters.

template<class Matrix, class Precond, class VectorRHS, class VectorX>
std::tuple<size_t, scalar_type> operator()(const Matrix &A, const Precond &P, const VectorRHS &rhs, const VectorX &x)

Computes the solution for the given system matrix A and the right-hand side rhs.

Returns the number of iterations made and the achieved relative residual as a std::tuple<size_t,scalar_type>. The solution vector x provides initial approximation on input and holds the computed solution on output.

template<class Precond, class VectorRHS, class VectorX>
std::tuple<size_t, scalar_type> operator()(const Precond &P, const VectorRHS &rhs, const VectorX &x)

Computes the solution for the given right-hand side rhs. The matrix that was used to create the preconditioner P is used as the system matrix.

Returns the number of iterations made and the achieved relative residual as a std::tuple<size_t,scalar_type>. The solution vector x provides initial approximation on input and holds the computed solution on output.

AMGCL implementats the following iterative solvers:

CG

template<class Backend, class InnerProduct = amgcl::detail::default_inner_product>
class amgcl::solver::cg

Include <amgcl/solver/cg.hpp>

The Conjugate Gradient method is an effective method for symmetric positive definite systems. It is probably the oldest and best known of the nonstationary methods [Barr94], [Saad03].

typedef typename Backend::value_type value_type

The value type of the system matrix

typedef typename amgcl::math::scalar_of<value_type>::type scalar_type

The scalar type corresponding to the value type. For example, when the value type is std::complex<double>, then the scalar type is double.

class params

The solver parameters.

size_t maxiter = 100

The maximum number of iterations

scalar_type tol = 1e-8

Target relative residual error \(\varepsilon = \frac{||f - Ax||}{|| f ||}\)

scalar_type abstol = std::numeric_limits<scalar_type>::min()

Target absolute residual error \(\varepsilon = ||f - Ax||\)

bool ns_search = false

Ignore the trivial solution x=0 when the RHS is zero. Useful when searching for the null-space vectors of the system.

bool verbose = false

Output the current iteration number and relative residual during solution.

BiCGStab

template<class Backend, class InnerProduct = amgcl::detail::default_inner_product>
class amgcl::solver::bicgstab

Include <amgcl/solver/bicgstab.hpp>

The BiConjugate Gradient Stabilized method (BiCGStab) was developed to solve nonsymmetric linear systems while avoiding the often irregular convergence patterns of the Conjugate Gradient Squared method [Barr94].

typedef typename Backend::value_type value_type

The value type of the system matrix

typedef typename amgcl::math::scalar_of<value_type>::type scalar_type

The scalar type corresponding to the value type. For example, when the value type is std::complex<double>, then the scalar type is double.

class params

The solver parameters.

bool check_after = false

Always do at least one iteration

amgcl::preconditioner::side::type pside = amgcl::preconditioner::side::right

Preconditioner kind (left/right)

size_t maxiter = 100

The maximum number of iterations

scalar_type tol = 1e-8

Target relative residual error \(\varepsilon = \frac{||f - Ax||}{|| f ||}\)

scalar_type abstol = std::numeric_limits<scalar_type>::min()

Target absolute residual error \(\varepsilon = ||f - Ax||\)

bool ns_search = false

Ignore the trivial solution x=0 when the RHS is zero. Useful when searching for the null-space vectors of the system.

bool verbose = false

Output the current iteration number and relative residual during solution.

BiCGStab(L)

template<class Backend, class InnerProduct = amgcl::detail::default_inner_product>
class amgcl::solver::bicgstabl

Include <amgcl/solver/bicgstabl.hpp>

This is a generalization of the BiCGStab method [SlDi93], [Fokk96]. For \(L=1\), this algorithm coincides with BiCGStab. In some situations it may be profitable to take \(L>2\). Although the steps of BiCGStab(L) are more expensive for larger L, numerical experiments indicate that, in certain situations, due to a faster convergence, for instance, BiCGStab(4) performs better than BiCGStab(2).

typedef typename Backend::value_type value_type

The value type of the system matrix

typedef typename amgcl::math::scalar_of<value_type>::type scalar_type

The scalar type corresponding to the value type. For example, when the value type is std::complex<double>, then the scalar type is double.

class params

The solver parameters.

int L = 2

The order of the method

scalar_type delta = 0

Threshold used to decide when to refresh computed residuals.

bool convex = true

Use a convex function of the MinRes and OR polynomials after the BiCG step instead of default MinRes

amgcl::preconditioner::side::type pside = amgcl::preconditioner::side::right

Preconditioner kind (left/right)

size_t maxiter = 100

The maximum number of iterations

scalar_type tol = 1e-8

Target relative residual error \(\varepsilon = \frac{||f - Ax||}{|| f ||}\)

scalar_type abstol = std::numeric_limits<scalar_type>::min()

Target absolute residual error \(\varepsilon = ||f - Ax||\)

bool ns_search = false

Ignore the trivial solution x=0 when the RHS is zero. Useful when searching for the null-space vectors of the system.

bool verbose = false

Output the current iteration number and relative residual during solution.

GMRES

template<class Backend, class InnerProduct = amgcl::detail::default_inner_product>
class amgcl::solver::gmres

Include <amgcl/solver/gmres.hpp>

The Generalized Minimal Residual method is an extension of MINRES (which is only applicable to symmetric systems) to unsymmetric systems. Like MINRES, it generates a sequence of orthogonal vectors, but in the absence of symmetry this can no longer be done with short recurrences; instead, all previously computed vectors in the orthogonal sequence have to be retained. For this reason, “restarted” versions of the method are used [Barr94].

typedef typename Backend::value_type value_type

The value type of the system matrix

typedef typename amgcl::math::scalar_of<value_type>::type scalar_type

The scalar type corresponding to the value type. For example, when the value type is std::complex<double>, then the scalar type is double.

class params

The solver parameters.

int M = 30

The number of iterations before restart

amgcl::preconditioner::side::type pside = amgcl::preconditioner::side::right

Preconditioner kind (left/right)

size_t maxiter = 100

The maximum number of iterations

scalar_type tol = 1e-8

Target relative residual error \(\varepsilon = \frac{||f - Ax||}{|| f ||}\)

scalar_type abstol = std::numeric_limits<scalar_type>::min()

Target absolute residual error \(\varepsilon = ||f - Ax||\)

bool ns_search = false

Ignore the trivial solution x=0 when the RHS is zero. Useful when searching for the null-space vectors of the system.

bool verbose = false

Output the current iteration number and relative residual during solution.

“Loose” GMRES (LGMRES)

template<class Backend, class InnerProduct = amgcl::detail::default_inner_product>
class amgcl::solver::lgmres

Include <amgcl/solver/lgmres.hpp>

The residual vectors at the end of each restart cycle of restarted GMRES often alternate direction in a cyclic fashion, thereby slowing convergence. LGMRES is an implementation of a technique for accelerating the convergence of restarted GMRES by disrupting this alternating pattern. The new algorithm resembles a full conjugate gradient method with polynomial preconditioning, and its implementation requires minimal changes to the standard restarted GMRES algorithm [BaJM05].

typedef typename Backend::value_type value_type

The value type of the system matrix

typedef typename amgcl::math::scalar_of<value_type>::type scalar_type

The scalar type corresponding to the value type. For example, when the value type is std::complex<double>, then the scalar type is double.

class params

The solver parameters.

unsigned K = 3

Number of vectors to carry between inner GMRES iterations. According to [BaJM05], good values are in the range of 1-3. However, if you want to use the additional vectors to accelerate solving multiple similar problems, larger values may be beneficial.

bool always_reset = true

Reset augmented vectors between solves. If the solver is used to repeatedly solve similar problems, then keeping the augmented vectors between solves may speed up subsequent solves. This flag, when set, resets the augmented vectors at the beginning of each solve.

int M = 30

The number of iterations before restart

amgcl::preconditioner::side::type pside = amgcl::preconditioner::side::right

Preconditioner kind (left/right)

size_t maxiter = 100

The maximum number of iterations

scalar_type tol = 1e-8

Target relative residual error \(\varepsilon = \frac{||f - Ax||}{|| f ||}\)

scalar_type abstol = std::numeric_limits<scalar_type>::min()

Target absolute residual error \(\varepsilon = ||f - Ax||\)

bool ns_search = false

Ignore the trivial solution x=0 when the RHS is zero. Useful when searching for the null-space vectors of the system.

bool verbose = false

Output the current iteration number and relative residual during solution.

Flexible GMRES (FGMRES)

template<class Backend, class InnerProduct = amgcl::detail::default_inner_product>
class amgcl::solver::fgmres

Include <amgcl/solver/fgmres.hpp>

Often, the application of the preconditioner P is a result of some unspecified computation, possibly another iterative process. In such cases, it may well happen that P is not a constant operator. The preconditioned iterative solvers may not converge if P is not constant. There are a number of variants of iterative procedures developed in the literature that can accommodate variations in the preconditioner, i.e., that allow the preconditioner to vary from step to step. Such iterative procedures are called “flexible” iterations. The method implements flexible variant of the GMRES algorithm [Saad03].

typedef typename Backend::value_type value_type

The value type of the system matrix

typedef typename amgcl::math::scalar_of<value_type>::type scalar_type

The scalar type corresponding to the value type. For example, when the value type is std::complex<double>, then the scalar type is double.

class params

The solver parameters.

int M = 30

The number of iterations before restart

size_t maxiter = 100

The maximum number of iterations

scalar_type tol = 1e-8

Target relative residual error \(\varepsilon = \frac{||f - Ax||}{|| f ||}\)

scalar_type abstol = std::numeric_limits<scalar_type>::min()

Target absolute residual error \(\varepsilon = ||f - Ax||\)

bool ns_search = false

Ignore the trivial solution x=0 when the RHS is zero. Useful when searching for the null-space vectors of the system.

bool verbose = false

Output the current iteration number and relative residual during solution.

IDR(s)

template<class Backend, class InnerProduct = amgcl::detail::default_inner_product>
class amgcl::solver::idrs

Include <amgcl/solver/idrs.hpp>

This is a very stable and efficient IDR(s) variant as described in [GiSo11]. The Induced Dimension Reduction method, IDR(s), is a robust and efficient short-recurrence Krylov subspace method for solving large nonsymmetric systems of linear equations.

IDR(s) compared to BI-CGSTAB/BiCGStab():

• Faster.
• More robust.
• More flexible.

typedef typename Backend::value_type value_type

The value type of the system matrix

typedef typename amgcl::math::scalar_of<value_type>::type scalar_type

The scalar type corresponding to the value type. For example, when the value type is std::complex<double>, then the scalar type is double.

class params

The solver parameters.

unsigned s = 4

Dimension of the shadow space in IDR(s).

scalar_type omega = 0.7

Computation of omega.

• If omega = 0, a standard minimum residual step is performed.
• If omega > 0, omega is increased if the cosine of the angle between Ar and r < omega.

bool smoothing = false

Specifies if residual smoothing must be applied.

bool replacement = false

Residual replacement. Determines the residual replacement strategy. If true, the recursively computed residual is replaced by the true residual.

size_t maxiter = 100

The maximum number of iterations

scalar_type tol = 1e-8

Target relative residual error \(\varepsilon = \frac{||f - Ax||}{|| f ||}\)

scalar_type abstol = std::numeric_limits<scalar_type>::min()

Target absolute residual error \(\varepsilon = ||f - Ax||\)

bool ns_search = false

Ignore the trivial solution x=0 when the RHS is zero. Useful when searching for the null-space vectors of the system.

bool verbose = false

Output the current iteration number and relative residual during solution.

Richardson iteration

template<class Backend, class InnerProduct = amgcl::detail::default_inner_product>
class amgcl::solver::richardson

Include <amgcl/solver/richardson.hpp>

The preconditioned Richardson iterative method

\[x^{i+1} = x^{i} + \omega P(f - A x^{i})\]

typedef typename Backend::value_type value_type

The value type of the system matrix

typedef typename amgcl::math::scalar_of<value_type>::type scalar_type

The scalar type corresponding to the value type. For example, when the value type is std::complex<double>, then the scalar type is double.

class params

The solver parameters.

scalar_type damping = 1.0

The damping factor \(\omega\)

size_t maxiter = 100

The maximum number of iterations

scalar_type tol = 1e-8

Target relative residual error \(\varepsilon = \frac{||f - Ax||}{|| f ||}\)

scalar_type abstol = std::numeric_limits<scalar_type>::min()

Target absolute residual error \(\varepsilon = ||f - Ax||\)

bool ns_search = false

Ignore the trivial solution x=0 when the RHS is zero. Useful when searching for the null-space vectors of the system.

bool verbose = false

Output the current iteration number and relative residual during solution.

PreOnly

template<class Backend, class InnerProduct = amgcl::detail::default_inner_product>
class amgcl::solver::preonly

Include <amgcl/solver/preonly.hpp>

Only apply the preconditioner once. This is not very useful as a standalone solver, but may be used in composite preconditioners as a nested solver, so that the composite preconditioner itself remains linear and may be used with a non-flexible iterative solver.

typedef typename Backend::value_type value_type

The value type of the system matrix

typedef typename amgcl::math::scalar_of<value_type>::type scalar_type

The scalar type corresponding to the value type. For example, when the value type is std::complex<double>, then the scalar type is double.

class params

The solver parameters.