Preconditioners

In this chapter, all preconditioners are presented. All preconditioners support local operators. They can be used as a global preconditioner via block-jacobi scheme which works locally on each interior matrix. To provide fast application, all preconditioners require extra memory to keep the approximated operator.

template<class OperatorType, class VectorType, typename ValueType>
class Preconditioner : public rocalution::Solver<OperatorType, VectorType, ValueType>

Base class for all preconditioners.

Template Parameters:

• OperatorType – - can be LocalMatrix or GlobalMatrix

• VectorType – - can be LocalVector or GlobalVector

• ValueType – - can be float, double, std::complex<float> or std::complex<double>

Subclassed by rocalution::AIChebyshev< OperatorType, VectorType, ValueType >, rocalution::AS< OperatorType, VectorType, ValueType >, rocalution::BlockJacobi< OperatorType, VectorType, ValueType >, rocalution::BlockPreconditioner< OperatorType, VectorType, ValueType >, rocalution::DiagJacobiSaddlePointPrecond< OperatorType, VectorType, ValueType >, rocalution::FSAI< OperatorType, VectorType, ValueType >, rocalution::GS< OperatorType, VectorType, ValueType >, rocalution::IC< OperatorType, VectorType, ValueType >, rocalution::ILU< OperatorType, VectorType, ValueType >, rocalution::ILUT< OperatorType, VectorType, ValueType >, rocalution::Jacobi< OperatorType, VectorType, ValueType >, rocalution::MultiColored< OperatorType, VectorType, ValueType >, rocalution::MultiElimination< OperatorType, VectorType, ValueType >, rocalution::SGS< OperatorType, VectorType, ValueType >, rocalution::SPAI< OperatorType, VectorType, ValueType >, rocalution::TNS< OperatorType, VectorType, ValueType >, rocalution::VariablePreconditioner< OperatorType, VectorType, ValueType >

Code Structure

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

Jacobi Method

template<class OperatorType, class VectorType, typename ValueType>
class Jacobi : public rocalution::Preconditioner<OperatorType, VectorType, ValueType>

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) \]

Template Parameters:

• OperatorType – - can be LocalMatrix or GlobalMatrix

• VectorType – - can be LocalVector or GlobalVector

• ValueType – - can be float, double, std::complex<float> or std::complex<double>

Note

Damping parameter \(\omega\) can be adjusted by rocalution::FixedPoint::SetRelaxation().

(Symmetric) Gauss-Seidel / (S)SOR Method

template<class OperatorType, class VectorType, typename ValueType>
class GS : public rocalution::Preconditioner<OperatorType, VectorType, ValueType>

Gauss-Seidel / Successive Over-Relaxation 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)\).

Template Parameters:

• OperatorType – - can be LocalMatrix

• VectorType – - can be LocalVector

• ValueType – - can be float, double, std::complex<float> or std::complex<double>

template<class OperatorType, class VectorType, typename ValueType>
class SGS : public rocalution::Preconditioner<OperatorType, VectorType, ValueType>

Symmetric Gauss-Seidel / Symmetric Successive Over-Relaxation Method.

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

Template Parameters:

• OperatorType – - can be LocalMatrix

• VectorType – - can be LocalVector

• ValueType – - can be float, double, std::complex<float> or std::complex<double>

Note

Relaxation parameter \(\omega\) can be adjusted by rocalution::FixedPoint::SetRelaxation().

Incomplete Factorizations

ILU

template<class OperatorType, class VectorType, typename ValueType>
class ILU : public rocalution::Preconditioner<OperatorType, VectorType, ValueType>

Incomplete LU Factorization based on levels.

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

Template Parameters:

• OperatorType – - can be LocalMatrix

• VectorType – - can be LocalVector

• ValueType – - can be float, double, std::complex<float> or std::complex<double>

virtual void rocalution::ILU::Set(int p, bool level = true)

Initialize ILU(p) factorization.

Initialize ILU(p) factorization based on power. [11]

• level = true build the structure based on levels

• level = false build the structure only based on the power(p+1)

ILUT

template<class OperatorType, class VectorType, typename ValueType>
class ILUT : public rocalution::Preconditioner<OperatorType, VectorType, ValueType>

Incomplete LU Factorization based on threshold.

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. [11]

Template Parameters:

• OperatorType – - can be LocalMatrix

• VectorType – - can be LocalVector

• ValueType – - can be float, double, std::complex<float> or std::complex<double>

virtual void rocalution::ILUT::Set(double t)

Set drop-off threshold.

virtual void rocalution::ILUT::Set(double t, int maxrow)

Set drop-off threshold and maximum fill-ins per row.

IC

template<class OperatorType, class VectorType, typename ValueType>
class IC : public rocalution::Preconditioner<OperatorType, VectorType, ValueType>

Incomplete Cholesky Factorization without fill-ins.

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.

Template Parameters:

• OperatorType – - can be LocalMatrix

• VectorType – - can be LocalVector

• ValueType – - can be float, double, std::complex<float> or std::complex<double>

AI Chebyshev

template<class OperatorType, class VectorType, typename ValueType>
class AIChebyshev : public rocalution::Preconditioner<OperatorType, VectorType, ValueType>

Approximate Inverse - Chebyshev Preconditioner.

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

Template Parameters:

• OperatorType – - can be LocalMatrix

• VectorType – - can be LocalVector

• ValueType – - can be float, double, std::complex<float> or std::complex<double>

void rocalution::AIChebyshev::Set(int p, ValueType lambda_min, ValueType lambda_max)

Set order, min and max eigenvalues.

FSAI

template<class OperatorType, class VectorType, typename ValueType>
class FSAI : public rocalution::Preconditioner<OperatorType, VectorType, ValueType>

Factorized Approximate Inverse Preconditioner.

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. [6]

Note

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

Template Parameters:

• OperatorType – - can be LocalMatrix

• VectorType – - can be LocalVector

• ValueType – - can be float, double, std::complex<float> or std::complex<double>

void rocalution::FSAI::Set(int power)

Set the power of the system matrix sparsity pattern.

void rocalution::FSAI::Set(const OperatorType &pattern)

Set an external sparsity pattern.

void rocalution::FSAI::SetPrecondMatrixFormat(unsigned int mat_format, int blockdim = 1)

Set the matrix format of the preconditioner.

SPAI

template<class OperatorType, class VectorType, typename ValueType>
class SPAI : public rocalution::Preconditioner<OperatorType, VectorType, ValueType>

SParse Approximate Inverse Preconditioner.

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}\). [5]

Template Parameters:

• OperatorType – - can be LocalMatrix

• VectorType – - can be LocalVector

• ValueType – - can be float, double, std::complex<float> or std::complex<double>

void rocalution::SPAI::SetPrecondMatrixFormat(unsigned int mat_format, int blockdim = 1)

Set the matrix format of the preconditioner.

TNS

template<class OperatorType, class VectorType, typename ValueType>
class TNS : public rocalution::Preconditioner<OperatorType, VectorType, ValueType>

Truncated Neumann Series Preconditioner.

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 from, the application of the preconditioner is based on matrix-vector operations only. The matrix format for the stored matrices can be specified.

Template Parameters:

• OperatorType – - can be LocalMatrix

• VectorType – - can be LocalVector

• ValueType – - can be float, double, std::complex<float> or std::complex<double>

void rocalution::TNS::Set(bool imp)

Set implicit (true) or explicit (false) computation.

void rocalution::TNS::SetPrecondMatrixFormat(unsigned int mat_format, int blockdim = 1)

Set the matrix format of the preconditioner.

MultiColored Preconditioners

template<class OperatorType, class VectorType, typename ValueType>
class MultiColored : public rocalution::Preconditioner<OperatorType, VectorType, ValueType>

Base class for all multi-colored preconditioners.

Template Parameters:

• OperatorType – - can be LocalMatrix

• VectorType – - can be LocalVector

• ValueType – - can be float, double, std::complex<float> or std::complex<double>

Subclassed by rocalution::MultiColoredILU< OperatorType, VectorType, ValueType >, rocalution::MultiColoredSGS< OperatorType, VectorType, ValueType >

void rocalution::MultiColored::SetPrecondMatrixFormat(unsigned int mat_format, int blockdim = 1)

Set a specific matrix type of the decomposed block matrices.

void rocalution::MultiColored::SetDecomposition(bool decomp)

Set if the preconditioner should be decomposed or not.

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

template<class OperatorType, class VectorType, typename ValueType>
class MultiColoredGS : public rocalution::MultiColoredSGS<OperatorType, VectorType, ValueType>

Multi-Colored Gauss-Seidel / SOR Preconditioner.

The Multi-Colored Symmetric 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 GS preconditioner.

Template Parameters:

• OperatorType – - can be LocalMatrix

• VectorType – - can be LocalVector

• ValueType – - can be float, double, std::complex<float> or std::complex<double>

template<class OperatorType, class VectorType, typename ValueType>
class MultiColoredSGS : public rocalution::MultiColored<OperatorType, VectorType, ValueType>

Multi-Colored Symmetric Gauss-Seidel / SSOR Preconditioner.

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 SGS preconditioner.

Template Parameters:

• OperatorType – - can be LocalMatrix

• VectorType – - can be LocalVector

• ValueType – - can be float, double, std::complex<float> or std::complex<double>

Subclassed by rocalution::MultiColoredGS< OperatorType, VectorType, ValueType >

void rocalution::MultiColoredSGS::SetRelaxation(ValueType omega)

Set the relaxation parameter for the SOR/SSOR scheme.

Note

The preconditioner matrix format can be changed using rocalution::MultiColored::SetPrecondMatrixFormat().

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

template<class OperatorType, class VectorType, typename ValueType>
class MultiColoredILU : public rocalution::MultiColored<OperatorType, VectorType, ValueType>

Multi-Colored Incomplete LU Factorization Preconditioner.

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. [8]

Template Parameters:

• OperatorType – - can be LocalMatrix

• VectorType – - can be LocalVector

• ValueType – - can be float, double, std::complex<float> or std::complex<double>

void rocalution::MultiColoredILU::Set(int p)

Initialize a multi-colored ILU(p, p+1) preconditioner.

void rocalution::MultiColoredILU::Set(int p, int q, bool level = true)

Initialize a multi-colored ILU(p, q) preconditioner.

level = true will perform the factorization with levels

level = false will perform the factorization only on the power(q)-pattern

Note

The preconditioner matrix format can be changed using rocalution::MultiColored::SetPrecondMatrixFormat().

Multi-Elimination Incomplete LU

template<class OperatorType, class VectorType, typename ValueType>
class MultiElimination : public rocalution::Preconditioner<OperatorType, VectorType, ValueType>

Multi-Elimination Incomplete LU Factorization Preconditioner.

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. [11]

Template Parameters:

• OperatorType – - can be LocalMatrix

• VectorType – - can be LocalVector

• ValueType – - can be float, double, std::complex<float> or std::complex<double>

inline int rocalution::MultiElimination::GetSizeDiagBlock(void) const

Returns the size of the first (diagonal) block of the preconditioner.

inline int rocalution::MultiElimination::GetLevel(void) const

Return the depth of the current level.

void rocalution::MultiElimination::Set(Solver<OperatorType, VectorType, ValueType> &AA_Solver, int level, double drop_off = 0.0)

Initialize (recursively) ME-ILU with level (depth of recursion)

AA_Solvers - defines the last-block solver

drop_off - defines drop-off tolerance

void rocalution::MultiElimination::SetPrecondMatrixFormat(unsigned int mat_format, int blockdim = 1)

Set a specific matrix type of the decomposed block matrices.

Diagonal Preconditioner for Saddle-Point Problems

template<class OperatorType, class VectorType, typename ValueType>
class DiagJacobiSaddlePointPrecond : public rocalution::Preconditioner<OperatorType, VectorType, ValueType>

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\).

Template Parameters:

• OperatorType – - can be LocalMatrix

• VectorType – - can be LocalVector

• ValueType – - can be float, double, std::complex<float> or std::complex<double>

void rocalution::DiagJacobiSaddlePointPrecond::Set(Solver<OperatorType, VectorType, ValueType> &K_Solver, Solver<OperatorType, VectorType, ValueType> &S_Solver)

Initialize solver for \(K\) and \(S\).

(Restricted) Additive Schwarz Preconditioner

template<class OperatorType, class VectorType, typename ValueType>
class AS : public rocalution::Preconditioner<OperatorType, VectorType, ValueType>

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 wich overlap. This leads to contributions from two preconditioners on the overlapped area which are scaled by \(1/2\). [2]

Template Parameters:

• OperatorType – - can be LocalMatrix

• VectorType – - can be LocalVector

• ValueType – - can be float, double, std::complex<float> or std::complex<double>

Subclassed by rocalution::RAS< OperatorType, VectorType, ValueType >

void rocalution::AS::Set(int nb, int overlap, Solver<OperatorType, VectorType, ValueType> **preconds)

Set number of blocks, overlap and array of preconditioners.

template<class OperatorType, class VectorType, typename ValueType>
class RAS : public rocalution::AS<OperatorType, VectorType, ValueType>

Restricted Additive Schwarz Preconditioner.

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. [2]

Template Parameters:

• OperatorType – - can be LocalMatrix

• VectorType – - can be LocalVector

• ValueType – - can be float, double, std::complex<float> or std::complex<double>

The overlapped area is shown in Fig. 13.

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

Block-Jacobi (MPI) Preconditioner

template<class OperatorType, class VectorType, typename ValueType>
class BlockJacobi : public rocalution::Preconditioner<OperatorType, VectorType, ValueType>

Block-Jacobi 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.

Template Parameters:

• OperatorType – - can be GlobalMatrix

• VectorType – - can be GlobalVector

• ValueType – - can be float, double, std::complex<float> or std::complex<double>

void rocalution::BlockJacobi::Set(Solver<LocalMatrix<ValueType>, LocalVector<ValueType>, ValueType> &precond)

Set local preconditioner.

The Block-Jacobi (MPI) preconditioner is shown in Fig. 14.

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

Block Preconditioner

template<class OperatorType, class VectorType, typename ValueType>
class BlockPreconditioner : public rocalution::Preconditioner<OperatorType, VectorType, ValueType>

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 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 SetDiagonalSolver().

Template Parameters:

• OperatorType – - can be LocalMatrix

• VectorType – - can be LocalVector

• ValueType – - can be float, double, std::complex<float> or std::complex<double>

void rocalution::BlockPreconditioner::Set(int n, const int *size, Solver<OperatorType, VectorType, ValueType> **D_solver)

Set number, size and diagonal solver.

void rocalution::BlockPreconditioner::SetDiagonalSolver(void)

Set diagonal solver mode.

void rocalution::BlockPreconditioner::SetLSolver(void)

Set lower triangular sweep mode.

void rocalution::BlockPreconditioner::SetExternalLastMatrix(const OperatorType &mat)

Set external last block matrix.

virtual void rocalution::BlockPreconditioner::SetPermutation(const LocalVector<int> &perm)

Set permutation vector.

Variable Preconditioner

template<class OperatorType, class VectorType, typename ValueType>
class VariablePreconditioner : public rocalution::Preconditioner<OperatorType, VectorType, ValueType>

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.

Template Parameters:

• OperatorType – - can be LocalMatrix

• VectorType – - can be LocalVector

• ValueType – - can be float, double, std::complex<float> or std::complex<double>

virtual void rocalution::VariablePreconditioner::SetPreconditioner(int n, Solver<OperatorType, VectorType, ValueType> **precond)

Set the preconditioner sequence.