@@ -113,7 +113,7 @@ The computational parameters that varied in the benchmark are:
\end{itemize}
To ensure comparable results obtained with different preconditioners for the BiCGstab method, care must be taken to select suitable stopping criteria for the iterative algorithm.
The linear system solvers implemented in TNL use left-preconditioning, whereas the Hypre library implements solvers with right-preconditioning (see \cref{chapter:linear systems}\todo{make sure left- vs right-preconditioning is described there} for details).
The linear system solvers implemented in TNL use left-preconditioning, whereas the Hypre library implements solvers with right-preconditioning (see \cref{chapter:linear systems} for details).
Recall that right-preconditioned methods operate with the \emph{original} residual and right-hand-side vectors, whereas left-preconditioned methods operate with \emph{preconditioned} residual and right-hand-side vectors.
Hence, assuming the Jacobi preconditioner, the right-hand-side vector is scaled with the matrix diagonal entries during left-preconditioning and the stopping criterion based on the norm of the residual vector divided by the norm of the right-hand-side vector cannot be used with the same threshold as for the unpreconditioned or right-preconditioned system.
To remedy this problem, we assemble the linear system with scaled rows to obtain ones on the main diagonal of the matrix (i.e., the Jacobi preconditioning is applied during the matrix assembly).
@@ -33,16 +33,86 @@ The induced dimension reduction methods, denoted as IDR($s$), are a new family o
Although they were reported to outperform traditional methods such as BiCGstab for several problems, they are still much less popular.
In general, the research focuses on the development of more efficient preconditioners rather than new iterative methods themselves.
\inline{domain decomposition and multigrid methods}
\inline{stopping criteria?}
\section{Preconditioning techniques}
\label{sec:preconditioning}
\inline{general intro, discuss left-preconditioning vs right-preconditioning}
\inline{make a summary of the state of the art}
\inline{stationary methods (Jacobi, SOR, SSOR), incomplete factorizations (IC and ILU), polynomial, sparse approximate inverse, multigrid (MG, problem-specific) or algebraic multigrid (AMG), GenEO, ...}
Preconditioning is the main technique to improve the robustness and efficiency of iterative methods.
The technique transforms the original linear system into a different system with the same solution, but which may be easier to solve using an iterative method.
Given a generic preconditioner represented by a non-singular matrix $\matrix M$, the original linear system
\begin{equation}
\matrix A \vec x = \vec b
\end{equation}
can be preconditioned either from the left, leading to the preconditioned system
\begin{equation}
\matrix M^{-1}\matrix A \vec x = \matrix M^{-1}\vec b,
\end{equation}
or from the right, leading to
\begin{equation}
\matrix A \matrix M^{-1}\vec y = \vec b, \qquad\vec x = \matrix M^{-1}\vec y,
\end{equation}
where the system is solved in terms of the transformed unknown variables $\vec y =\matrix M \vec x$.
Finally, a split-preconditioning in the form
\begin{equation}
\matrix M_1^{-1}\matrix A \matrix M_2^{-1}\vec y = \matrix M_1^{-1}\vec b, \qquad\vec x = \matrix M_2^{-1}\vec y
\end{equation}
is also possible, where the preconditioner is $\matrix M =\matrix M_1\matrix M_2$.
There is no consensus on which kind of preconditioning is the best and how to choose the preconditioning matrix $\matrix M$.
The choices typically depend on the iterative method used as well as on specific matrix properties for the given problem, however, several general-purpose preconditioners that do not require problem-specific input have been developed.
In general, the minimal requirements imposed on the preconditioner $\matrix M$ are that it must be non-singular and it should approximate $\matrix A$ in some sense, such that $\matrix M^{-1}\matrix A$ and $\matrix A \matrix M^{-1}$ are close to the identity matrix $\matrix I$.
Note that the preconditioned matrices $\matrix M^{-1}\matrix A$ and $\matrix A \matrix M^{-1}$ are not assembled explicitly as that would be too expensive and sparsity would be lost.
Instead, additional requirements are that the operation $\matrix M^{-1}\vec v$ should be easy to apply to an arbitrary vector $\vec v$ and it should offer opportunities for an efficient parallel implementation.
However, in practice, the last three requirements often go against each other, which poses a major difficulty to the development of high performance preconditioners.
\subsection{State of the art}
There are two general approaches to constructing preconditioners.
One approach is to design specialized algorithms that are (nearly) optimal for a small class of problems.
These typically utilize the knowledge of the underlying problem, such as the original partial differential equation from which the discrete problem arose.
While these preconditioners are often very effective, obtaining or using the underlying information may not be feasible or desirable.
Hence, the other approach comprising general purpose, purely algebraic preconditioning techniques may be preferable.
These preconditioners, while not optimal for any specific problem, achieve reasonable efficiency for a large class of problems and they are often easier to implement and use in software frameworks.
In the following, we give a brief overview of the most successful techniques.
Detailed reviews are available in, e.g., \cite{benzi:2002,gambolati:2011,pearson:2020}.
The simplest preconditioning techniques are based on stationary methods such as Jacobi, Gauss--Seidel, SOR, or their symmetric variants \cite{saad:2003iterative,benzi:2002}.
Another traditional class of generic preconditioners are incomplete factorization techniques, such as incomplete Cholesky (IC) and incomplete LU (ILU), which bring ideas developed for direct solution of sparse linear systems to iterative methods.
There are multiple strategies to control the fill-in by selecting the sparsity pattern of the factors either statically based on the matrix entry positions or dynamically based on the matrix entry values \cite{saad:1994ilut,li:2003crout,bollhofer:2001,mittal:2003}.
Additionally, algebraic multilevel techniques originating from domain decomposition and multigrid methods were applied also to incomplete factorization methods \cite{saad:1999bilutm,botta:1999,henon:2006} and algorithmic variants were developed to address parallelization and scalability issues \cite{karypis:1997,hysom:1999,hysom:2001,wang:2003,anzt:2018parilut,coleman:2018}.
Several incomplete factorization algorithms were also ported to GPU accelerators \cite{naumov:2011,naumov:2015}, though their efficiency remains limited compared to other preconditioning techniques.
Sparse approximate inverse preconditioners provide an alternative approach to remedy problems arising when incomplete factorizations are applied to indefinite systems or systems that are not diagonally dominant \cite{benzi:2002,saad:2003iterative}.
They also attracted attention thanks to their potential in parallel environments \cite{tang:1999}, including GPU accelerators \cite{dehnavi:2013,oyarzun:2014,bertaccini:2016,bernaschi:2019}.
Approximate inverse preconditioners may be applied to symmetric as well as non-symmetric systems and they may be computed in a non-factored form where the preconditioner is expressed as a single matrix, or in a factored form where the preconditioner is expressed as a product of two or more matrices.
For both classes there exist multiple completely different approaches to compute the approximate inverse or approximate inverse factors, the two main approaches are Frobenius norm minimization \cite{chow:1994,chow:2000,kolotilina:1993,kolotilina:1995,kolotilina:2000,kolotilina:1999} and incomplete bi-conjugation \cite{benzi:1996,benzi:1998,bridson:2000,bertaccini:2016}.
The sparse approximate inverse algorithms based on bi-conjugation show algebraic behavior that is similar to the incomplete LU factorizations \cite{bollhofer:2002a,bollhofer:2002b}.
When matrix values are changed but the sparsity pattern remains the same, it is possible to reuse an existing approximate inverse preconditioner to compute the updated approximate inverse more efficiently \cite{kolotilina:2000,bridson:2000}.
Several dynamic pattern selection strategies and multilevel techniques for sparse approximate inverse preconditioners have been recently investigated \cite{janna:2011,janna:2015,bernaschi:2019,franceschini:2018}.
Polynomial preconditioning is closely related to the development of Krylov methods.
It utilizes spectral information of the linear system matrix and is favorable for parallel architectures where the sparse matrix--vector multiplication delivers high performance.
Several types of polynomials can be used, the Chebyshev \cite{ashby:1992,asgasri:2009}, least-squares \cite{vangijzen:1995,li:2013}, and Newton \cite{liu:2015} polynomials are the most common in practice.
Several non-algebraic techniques were developed for solving linear systems arising from the discretization of partial differential equations.
The convergence rate of Krylov subspace methods preconditioned by the aforementioned techniques applied to these linear systems deteriorates as the systems become larger \cite{greenbaum:1997iterative,saad:2003iterative}.
Consequently, the efficiency of these methods drops significantly due to increasing number of iterations combined with the sheer problem size.
To this effect, multigrid methods \cite{wesseling:1991,briggs:2000,trottenberg:2001} were developed to obtain (nearly) optimal techniques in the sense that the number of iterations does not depend on the problem size, or increases only mildly.
Although multigrid, hereafter referred to as geometric multigrid (GMG), was originally developed as a method for solving partial differential equations, it was later generalized into algebraic multigrid (AMG) techniques \cite{ruge-stuben:1987,vanek:1996,stuben:2001,trottenberg:2001,lottes:2017,xu:2017} that can be applied even to problems not involving partial differential equations at all.
Unfortunately, multigrid techniques are naturally communication-bound due to coarse grid computations, interpolation and restriction operations.
Consequently, both geometric multigrid \cite{bauer:2016,goddeke:2010,feng:2008,zheng:2014} and algebraic multigrid \cite{baker:2011,bell:2012,gandham:2014,liu_yang_chen:2015} have been recently subject of intensive research to address the efficiency and scalability challenges on modern parallel platforms, including GPU accelerators.
Domain decomposition methods \cite{smith:1996,toselli:2004} are another class of techniques developed to efficiently solve partial differential equations.
Contrary to multigrid, domain decomposition methods were associated with parallel processing since their origin and massively parallel implementations for contemporary platforms exist of both main domain decomposition methods, BDDC and FETI \cite{sousedik:2013,kozubek:2013}.
Although both multigrid and domain decomposition were derived as methods for the solution of partial differential equations, they are nowadays used mostly as preconditioners to accelerate the convergence of Krylov subspace methods.
Multigrid and domain decomposition methods borrowed ideas from each other during their development, which lead to the development of various algebraic multilevel preconditioning techniques \cite{axelsson:1989,axelsson:1990,dryja:1996}.
Similarly to multigrid, recent development of domain decomposition methods tends towards purely algebraic preconditioning techniques \cite{agullo:2013,tselepidis:2019}.
To conclude the section, we mention some general techniques that can be used when multiple preconditioners accelerating the solution of a linear system are available.
Rather than formulating a new preconditioner that combines the features of multiple preconditioners, they can be applied together within a \emph{multi-preconditioned} iterative method \cite{bridson:2006,spillane:2016,greif:2017}.
More generally, when solving sequences of similar linear systems, such as those arising from the discretization of transient partial differential equations, it may be possible to reuse information from previously computed problems, including preconditioners \cite{simoncini:2007,carr2021recycling,amritkar2015recycling,gaul2014recycling}.
