Relaxation

A relaxation method or a smoother is used on each level of the AMG hierarchy during solution phase.

Damped Jacobi

template<class Backend>
class amgcl::relaxation::damped_jacobi

Include <amgcl/relaxation/damped_jacobi.hpp>

The damped Jacobi relaxation

class params

Damped Jacobi relaxation parameters

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.

scalar_type damping = 0.72

The damping factor

Gauss-Seidel

template<class Backend>
class amgcl::relaxation::gauss_seidel

Include <amgcl/relaxation/gauss_seidel.hpp>

The Gauss-Seidel relaxation. The relaxation is only available for the backends where the matrix supports row-wise iteration over its non-zero values.

class params

Gauss-Seidel relaxation parameters

bool serial = false

Use the serial version of the algorithm

Chebyshev

template<class Backend>
class amgcl::relaxation::chebyshev

Include <amgcl/relaxation/chebyshev.hpp>

Chebyshev iteration is an iterative method for the solution of a system of linear equations, and unlike Jacobi, it is not a stationary method. However, it does not require inner products like many other nonstationary methods (most Krylov methods). These inner products can be a performance bottleneck on certain distributed memory architectures. Furthermore, Chebyshev iteration is, like Jacobi, easier to parallelize than for instance Gauss–Seidel smoothers. The Chebyshev iteration requires some information about the spectrum of the matrix. For symmetric matrices, it needs an upper bound for the largest eigenvalue and a lower bound for the smallest eigenvalue [GhKK12].

class params

Chebyshev relaxation parameters

unsigned degree = 5

The degree of Chebyshev polynomial

float higher = 1.0

The highest eigen value safety upscaling. use boosting factor for a more conservative upper bound estimate [ABHT03].

float lower = 1.0 / 30

Lowest-to-highest eigen value ratio.

int power_iters = 0

The number of power iterations to apply for the spectral radius estimation. When 0, use Gershgorin disk theorem to estimate the spectral radius.

bool scale = false

Scale the system matrix

Incomplete LU relaxation

The incomplete LU factorization process computes a sparse lower triangular matrix \(L\) and a sparse upper triangular matrix \(U\) so that the residual matrix \(R = LU - A\) satisfies certain constraints, such as having zero entries in some locations. The relaxations in this section use various approaches to computation of the triangular factors \(L\) and \(U\), but share the triangular system solution implementation required in order to apply the relaxation. The parameters for the triangular solution algorithm are defined as follows:

template<class Backend>
class amgcl::relaxation::detail::ilu_solve

For the builtin OpenMP backend the incomplete triangular factors are solved using the OpenMP-parallel level scheduling approach. For the GPGPU backends, the triangular factors are solved approximately, using multiple damped Jacobi iterations [ChPa15].

class params
bool serial = false

Use the serial version of the algorithm. This parameter is only used with the amgcl::backend::builtin backend.

unsigned iters = 2

The number of Jacobi iterations to approximate the triangular system solution. This parameter is only used with GPGPU backends.

scalar_type damping = 1.0

The damping factor for the triangular solve approximation. This parameter is only used with GPGPU backends.

ILU0

template<class Backend>
class amgcl::relaxation::ilu0

Include <amgcl/relaxation/ilu0.hpp>

The incomplete LU factorization with zero fill-in [Saad03]. The zero pattern for the triangular factors \(L\) and \(U\) is taken to be exactly the zero pattern of the system matrix \(A\).

class params

ILU0 relaxation parameters

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.

scalar_type damping = 1.0

The damping factor

typename amgcl::relaxation::detail::ilu_solve<Backend>::params solve

The parameters for the triangular factor solver

ILUK

template<class Backend>
class amgcl::relaxation::iluk

Include <amgcl/relaxation/iluk.hpp>

The ILU(k) relaxation.

The incomplete LU factorization with the level of fill-in [Saad03]. The accuracy of the ILU0 incomplete factorization may be insufficient to yield an adequate rate of convergence. More accurate incomplete LU factorizations are often more efficient as well as more reliable. These more accurate factorizations will differ from ILU(0) by allowing some fill-in. Thus, ILUK(k) keeps the ‘k-th order fill-ins’ [Saad03].

The ILU(1) factorization results from taking the zero pattern for triangular factors to be the zero pattern of the product \(L_0 U_0\) of the factors \(L_0\), \(U_0\) obtained from ILU(0). This process is repeated to obtain the higher level of fill-in factorizations.

class params

ILUK relaxation parameters

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.

int k = 1

The level of fill-in

scalar_type damping = 1.0

The damping factor

typename amgcl::relaxation::detail::ilu_solve<Backend>::params solve

The parameters for the triangular factor solver

ILUP

template<class Backend>
class amgcl::relaxation::ilup

Include <amgcl/relaxation/ilup.hpp>

The ILUP(k) relaxation.

This variant of the ILU relaxation is similar to ILUK, but differs in the way the zero pattern for the triangular factors is determined. Instead of the recursive definition using the product \(LU\) of the factors from the previous level of fill-in, ILUP uses the powers of the boolean matrix \(S\) sharing the zero pattern with the system matrix \(A\) [MiKu03]. ILUP(0) coinsides with ILU0, ILUP(1) has the same zero pattern as \(S^2\), etc.

class params

ILUP relaxation parameters

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.

int k = 1

The level of fill-in

scalar_type damping = 1.0

The damping factor

typename amgcl::relaxation::detail::ilu_solve<Backend>::params solve

The parameters for the triangular factor solver

ILUT

template<class Backend>
class amgcl::relaxation::ilut

Include <amgcl/relaxation/ilut.hpp>

The \(\mathrm{ILUT}(p,\tau)\) relaxation.

Incomplete factorizations which rely on the levels of fill are blind to numerical values because elements that are dropped depend only on the structure of A. This can cause some difficulties for realistic problems that arise in many applications. A few alternative methods are available which are based on dropping elements in the Gaussian elimination process according to their magnitude rather than their locations. With these techniques, the zero pattern P is determined dynamically.

A generic ILU algorithm with threshold can be derived from the IKJ version of Gaussian elimination by including a set of rules for dropping small elements. In the factorization \(\mathrm{ILUT}(p,\tau)\), the following rule is used:

  1. an element is dropped (i.e., replaced by zero) if it is less than the relative tolerance \(\tau_i\) obtained by multiplying \(\tau\) by the original 2-norm of the i-th row.
  2. Only the \(p l_i\) largest elements are kept in the \(L\) part of the row and the \(p u_i\) largest elements in the \(U\) part of the row in addition to the diagonal element, which is always kept. \(l_i\) and \(u_i\) are the number of nonzero elements in the i-th row of the system matrix \(A\) below and above the diagonal.
class params

ILUT relaxation parameters

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.

scalar_type p = 2

The fill factor

scalar_type tau = 1e-2

The minimum magnitude of non-zero elements relative to the current row norm.

scalar_type damping = 1.0

The damping factor

typename amgcl::relaxation::detail::ilu_solve<Backend>::params solve

The parameters for the triangular factor solver

Sparse Approximate Inverse relaxation

Sparse approximate inverse (SPAI) smoothers based on the SPAI algorithm by Grote and Huckle [GrHu97]. The SPAI algorithm computes an approximate inverse \(M\) explicitly by minimizing \(I - MA\) in the Frobenius norm. Both the computation of \(M\) and its application as a smoother are inherently parallel. Since an effective sparsity pattern of \(M\) is in general unknown a priori, the computation cost can be greately reduced by choosing an a priori sparsity pattern for \(M\). For SPAI-0 and SPAI-1 the sparsity pattern of \(M\) is fixed: \(M\) is diagonal for SPAI-0, whereas for SPAI-1 the sparsity pattern of \(M\) is that of \(A\) [BrGr02].

SPAI0

template<class Backend>
class amgcl::relaxation::spai0

Include <amgcl/relaxation/spai0.hpp>

The SPAI-0 variant of the sparse approximate inverse smother [BrGr02].

class params

The SPAI-0 has no parameters

SPAI1

template<class Backend>
class amgcl::relaxation::spai1

Include <amgcl/relaxation/spai1.hpp>

The SPAI-1 variant of the sparse approximate inverse smother [BrGr02].

class params

The SPAI-1 has no parameters

Scalar to Block convertor

template<class BlockBackend, template<class> class Relax>
class amgcl::relaxation::as_block

Include <amgcl/relaxation/as_block.hpp>

Wrapper for the specified relaxation. Converts the input matrix from scalar to block format before constructing an amgcl smoother. See the Using near null-space vectors tutorial.

template <class Backend>
class type

The resulting relaxation class.