rocALUTION Preconditioners

This document describes the theory and usage of rocALUTION preconditioners. All preconditioners support local operators. They can be used as a global preconditioner via a block-Jacobi scheme, which works locally on each interior matrix. To provide fast application, all preconditioners require extra memory to keep the approximated operator.

Member documentation, configuration routines, and usage notes for each class are on the rocALUTION API library page.

Code structure

The preconditioners provide a solution to the system \(Mz = r\), where the solution \(z\) is either directly computed by the approximation scheme or iteratively obtained with \(z = 0\) initial guess.

rocalution::Preconditioner

Jacobi method

The Jacobi method is for solving a diagonally dominant system of linear equations \(Ax=b\). It solves for each diagonal element iteratively until convergence, such that

\[x_{i}^{(k+1)} = (1 - \omega)x_{i}^{(k)} + \frac{\omega}{a_{ii}} \left( b_{i} - \sum\limits_{j=1}^{i-1}{a_{ij}x_{j}^{(k)}} - \sum\limits_{j=i}^{n}{a_{ij}x_{j}^{(k)}} \right)\]

Note

To adjust the damping parameter \(\omega\), use rocalution::FixedPoint::SetRelaxation().

rocalution::Jacobi

(Symmetric) Gauss-Seidel or (S)SOR method

The Gauss-Seidel / SOR method is for solving system of linear equations \(Ax=b\). It approximates the solution iteratively with

\[x_{i}^{(k+1)} = (1 - \omega) x_{i}^{(k)} + \frac{\omega}{a_{ii}} \left( b_{i} - \sum\limits_{j=1}^{i-1}{a_{ij}x_{j}^{(k+1)}} - \sum\limits_{j=i}^{n}{a_{ij}x_{j}^{(k)}} \right),\]

with \(\omega \in (0,2)\).

The Symmetric Gauss-Seidel / SSOR method is for solving system of linear equations \(Ax=b\). It approximates the solution iteratively.

Note

To adjust the relaxation parameter \(\omega\), use rocalution::FixedPoint::SetRelaxation().

rocalution::GS, rocalution::SGS

Incomplete factorizations

ILU

The Incomplete LU Factorization based on levels computes a sparse lower and sparse upper triangular matrix such that \(A = LU - R\).

rocalution::ILU

ILUT

The Incomplete LU Factorization based on threshold computes a sparse lower and sparse upper triangular matrix such that \(A = LU - R\). Fill-in values are dropped depending on a threshold and number of maximal fill-ins per row.

rocalution::ILUT

IC

The Incomplete Cholesky Factorization computes a sparse lower triangular matrix such that \(A=LL^{T} - R\). Additional fill-ins are dropped and the sparsity pattern of the original matrix is preserved.

rocalution::IC

AI Chebyshev

The Approximate Inverse - Chebyshev Preconditioner is an inverse matrix preconditioner with values from a linear combination of matrix-valued Chebyshev polynomials.

rocalution::AIChebyshev

FSAI

The Factorized Sparse Approximate Inverse preconditioner computes a direct approximation of \(M^{-1}\) by minimizing the Frobenius norm \(||I - GL||_{F}\), where \(L\) denotes the exact lower triangular part of \(A\) and \(G:=M^{-1}\). The FSAI preconditioner is initialized by \(q\), based on the sparsity pattern of \(|A^{q}|\). However, it is also possible to supply external sparsity patterns in form of the LocalMatrix class.

Note

The FSAI preconditioner is only suited for symmetric positive definite matrices.

rocalution::FSAI

SPAI

The SParse Approximate Inverse algorithm is an explicitly computed preconditioner for general sparse linear systems. In its current implementation, only the sparsity pattern of the system matrix is supported. The SPAI computation is based on the minimization of the Frobenius norm \(||AM - I||_{F}\).

rocalution::SPAI

TNS

The Truncated Neumann Series (TNS) preconditioner is based on \(M^{-1} = K^{T} D^{-1} K\), where \(K=(I-LD^{-1}+(LD^{-1})^{2})\), with the diagonal \(D\) of \(A\) and the strictly lower triangular part \(L\) of \(A\). The preconditioner can be computed in two forms - explicitly and implicitly. In the explicit form, the full construction of \(M\) is performed via matrix-matrix operations, whereas in the implicit form, the application of the preconditioner is based on matrix-vector operations only. The matrix format for the stored matrices can be specified.

rocalution::TNS

MultiColored preconditioners

Multi-colored preconditioners reorder the unknowns to expose parallelism in forward and backward substitution. Derived multi-colored preconditioners can change the preconditioner matrix format with rocalution::MultiColored::SetPrecondMatrixFormat().

rocalution::MultiColored

MultiColored (symmetric) Gauss-Seidel / (S)SOR

The Multi-Colored Symmetric Gauss-Seidel / SSOR preconditioner is based on the splitting of the original matrix. Higher parallelism in solving the forward and backward substitution is obtained by performing a multi-colored decomposition. Details on the Symmetric Gauss-Seidel / SSOR algorithm can be found in the rocalution::SGS preconditioner.

The Multi-Colored Gauss-Seidel / SOR preconditioner is based on the splitting of the original matrix. Higher parallelism in solving the forward substitution is obtained by performing a multi-colored decomposition. Details on the Gauss-Seidel / SOR algorithm can be found in the rocalution::GS preconditioner.

Note

To change the preconditioner matrix format, use rocalution::MultiColored::SetPrecondMatrixFormat().

rocalution::MultiColoredGS, rocalution::MultiColoredSGS

MultiColored power(q)-pattern method ILU(p,q)

Multi-Colored Incomplete LU Factorization based on the ILU(p) factorization with a power(q)-pattern method. This method provides a higher degree of parallelism of forward and backward substitution compared to the standard ILU(p) preconditioner.

Note

To change the preconditioner matrix format, use rocalution::MultiColored::SetPrecondMatrixFormat().

rocalution::MultiColoredILU

Multi-elimination incomplete LU

The Multi-Elimination Incomplete LU preconditioner is based on the following decomposition

\[\begin{split}A = \begin{pmatrix} D & F \\ E & C \end{pmatrix} = \begin{pmatrix} I & 0 \\ ED^{-1} & I \end{pmatrix} \times \begin{pmatrix} D & F \\ 0 & \hat{A} \end{pmatrix},\end{split}\]

where \(\hat{A} = C - ED^{-1} F\). To make the inversion of \(D\) easier, we permute the preconditioning before the factorization with a permutation \(P\) to obtain only diagonal elements in \(D\). The permutation here is based on a maximal independent set. This procedure can be applied to the block matrix \(\hat{A}\), in this way we can perform the factorization recursively. In the last level of the recursion, we need to provide a solution procedure. By the design of the library, this can be any kind of solver.

rocalution::MultiElimination

Diagonal preconditioner for saddle-point problems

Consider the following saddle-point problem

\[\begin{split}A = \begin{pmatrix} K & F \\ E & 0 \end{pmatrix}.\end{split}\]

For such problems we can construct a diagonal Jacobi-type preconditioner of type

\[\begin{split}P = \begin{pmatrix} K & 0 \\ 0 & S \end{pmatrix},\end{split}\]

with \(S=ED^{-1}F\), where \(D\) are the diagonal elements of \(K\). The matrix \(S\) is fully constructed (via sparse matrix-matrix multiplication). The preconditioner needs to be initialized with two external solvers/preconditioners - one for the matrix \(K\) and one for the matrix \(S\).

rocalution::DiagJacobiSaddlePointPrecond

(Restricted) Additive Schwarz preconditioner

The Additive Schwarz preconditioner relies on a preconditioning technique, where the linear system \(Ax=b\) can be decomposed into small sub-problems based on \(A_{i} = R_{i}^{T}AR_{i}\), where \(R_{i}\) are restriction operators. Those restriction operators produce sub-matrices which overlap. This leads to contributions from two preconditioners on the overlapped area which are scaled by \(1/2\).

The Restricted Additive Schwarz preconditioner relies on a preconditioning technique, where the linear system \(Ax=b\) can be decomposed into small sub-problems based on \(A_{i} = R_{i}^{T}AR_{i}\), where \(R_{i}\) are restriction operators. The RAS method is a mixture of block Jacobi and the AS scheme. In this case, the sub-matrices contain overlapped areas from other blocks, too.

See the overlapped area in the figure below:

Fig. 13 Example of a 4 block-decomposed matrix - Additive Schwarz with overlapping preconditioner (left) and Restricted Additive Schwarz preconditioner (right).

rocalution::AS, rocalution::RAS

Block-Jacobi (MPI) preconditioner

The Block-Jacobi preconditioner is designed to wrap any local preconditioner and apply it in a global block fashion locally on each interior matrix.

See the Block-Jacobi (MPI) preconditioner in the figure below:

Fig. 14 Example of a 4 block-decomposed matrix - Block-Jacobi preconditioner.

rocalution::BlockJacobi

Block preconditioner

When handling vector fields, typically one can try to use different preconditioners and/or solvers for the different blocks. For such problems, the library provides a block-type preconditioner. This preconditioner builds the following block-type matrix

\[\begin{split}P = \begin{pmatrix} A_{d} & 0 & . & 0 \\ B_{1} & B_{d} & . & 0 \\ . & . & . & . \\ Z_{1} & Z_{2} & . & Z_{d} \end{pmatrix}\end{split}\]

The solution of \(P\) can be performed in two ways. It can be solved by block-lower-triangular sweeps with inversion of the blocks \(A_{d} \ldots Z_{d}\) and with a multiplication of the corresponding blocks. This is set by rocalution::BlockPreconditioner::SetLSolver() (which is the default solution scheme). Alternatively, it can be used only with an inverse of the diagonal \(A_{d} \ldots Z_{d}\) (Block-Jacobi type) by using rocalution::BlockPreconditioner::SetDiagonalSolver().

rocalution::BlockPreconditioner

Variable preconditioner

The Variable Preconditioner can hold a selection of preconditioners. Thus, any type of preconditioners can be combined. As example, the variable preconditioner can combine Jacobi, GS and ILU - then, the first iteration of the iterative solver will apply Jacobi, the second iteration will apply GS and the third iteration will apply ILU. After that, the solver will start again with Jacobi, GS, ILU.

rocalution::VariablePreconditioner