Preconditioning for Sparse Linear Systems at the Dawn of the 21st Century: History, Current Developments, and Future Perspectives

Abstract

Iterative methods are currently the solvers of choice for large sparse linear systems of equations. However, it is well known that the key factor for accelerating, or even allowing for, convergence is the preconditioner. The research on preconditioning techniques has characterized the last two decades. Nowadays, there are a number of different options to be considered when choosing the most appropriate preconditioner for the specific problem at hand. The present work provides an overview of the most popular algorithms available today, emphasizing the respective merits and limitations. The overview is restricted to algebraic preconditioners, that is, general-purpose algorithms requiring the knowledge of the system matrix only, independently of the specific problem it arises from. Along with the traditional distinction between incomplete factorizations and approximate inverses, the most recent developments are considered, including the scalable multigrid and parallel approaches which represent the current frontier of research. A separate section devoted to saddle-point problems, which arise in many different applications, closes the paper.

1. Historical Background

The term “preconditioning” generally refers to “the art of transforming a problem that appears intractable into another whose solution can be approximated rapidly” [1], while the “preconditioner” is the mathematical operator responsible for such a transformation. In the context of linear systems of equations, 𝐴𝐱=𝐛, (1.1) where 𝐴∈ℝ^𝑛×𝑛 and 𝐱,𝐛∈ℝ^𝑛, the preconditioner is a matrix which transforms (1.1) into an equivalent problem for which the convergence of an iterative method is much faster. Preconditioners are generally used when the matrix 𝐴 is large and sparse, as it typically, but not exclusively, arises from the numerical discretization of Partial Differential Equations (PDEs). It is not rare that the preconditioner itself is the operator which makes it possible to solve numerically the system (1.1), that otherwise would be intractable from a computational viewpoint.

The development of preconditioning techniques for large sparse linear systems is strictly connected to the history of iterative methods. As mentioned by Benzi [2], the term “preconditioning” used to identify the acceleration of the iterative solution of a linear system can be found for the first time in a 1968 paper by Evans [3], but the concept is pretty older and was already introduced by Cesari back in 1937 [4]. However, it is well known that the earliest tricks to solve numerically a linear system were already developed in the 19th century. As reported in the survey by Saad and Van der Vorst [5] and according to the work by Varga [6], Gauss set the pace for iterative methods in an 1823 letter, where the famous German mathematician describes a number of clever tricks to solve a singular 4×4 linear system arising from a geodetic problem. Most of the strategies suggested by Gauss compose what we currently know as the stationary Gauss-Seidel method, which was later formalized by Seidel in 1874. In his letter, Gauss jokes saying that the method he suggests is actually a pleasant entertainment, as one can think about many other things while carrying out the repetitive operations required to solve the systems. Actually, it is rather a surprise that such a method attracted a great interest in the era of automatic computing! In 1846 another famous algorithm was introduced by Jacobi to solve a sequence of 7×7 linear systems enforcing diagonal dominance by plane rotations. Other similar methods were later developed, such as the Richardson [7] and Cimmino [8] iteration, but no significant improvements in the field of numerical solution of linear systems can be reported until the 40s. This period is referred to as the prehistory of this subject by Benzi [2].

As it often occurred in the past, a dramatic and tragic historical event like World War II increased the research fundings for military reasons and helped accelerate the development and progress also in the field of numerical linear algebra. World War II brought a great development of the first automatic computing machines, which later led to the modern digital electronic era. With the availability of such new computing tools, the interest in the numerical solution of relatively large linear systems of equations suddenly grew. The more significant advances were obtained in the field of direct methods. Between the 40s and 70s the algorithms of choice for solving automatically a linear system are based on the Gaussian elimination, with the development of several important variants which helped greatly improve the native method. First, it was recognized that pivoting techniques are of paramount importance to stabilize the numerical computations of the triangular factors of the system matrix and reduce the effects of rounding errors, for example, [9]. Second, several reordering techniques were developed in order to limit the fill-in of the triangular factors and the number of operations required for their computation. The most famous bandwidth reducing algorithms, such as Cuthill-McKee and its reverse variant [10], nested dissection [11–13], and minimum degree [14–16], are still very popular today. Third, appropriate scaling strategies were introduced so as to balance the system matrix entries and help stabilize the triangular factor computation, for example, [17, 18].

In the meantime, iterative methods lived in a kind of niche, despite the publication in 1952 of the famous papers by Hestenes and Stiefel [19] and Lanczos [20] who discovered independently and almost simultaneously the Conjugate Gradient (CG) method for solving symmetric positive definite (SPD) linear systems. These papers did not go unnoticed, but at that time CG was interpreted as a direct method converging in 𝑚 steps, if 𝑚 is the number of distinct eigenvalues of the system matrix 𝐴. As it was soon realized [21], this property is lost in practice when working in finite arithmetics, so that convergence is actually achieved in a number of iterations possibly much larger than 𝑚. This is why CG was considered an attractive alternative to Gaussian elimination in a few lucky cases only and was not rated much more than just an elegant theory. The iterative method of choice of the period was the Successive Overrelaxation (SOR) algorithm, which was introduced independently by Young [22] in his Ph.D. thesis and Frankel [23] as an acceleration of the old Gauss-Seidel stationary iteration. The possibility of estimating theoretically the optimal overrelaxation parameter for a large class of matrices having the so-called property A [22] gained a remarkable success for SOR especially in problems arising from the discretization of PDEs of elliptic type, for example, in groundwater hydrology or nuclear reactor diffusion [6].

In the 60s the Finite Element (FE) method was introduced in structural mechanics. The new method had a great success and gave rise to a novel set of large matrices which were very sparse but neither diagonally dominant nor characterized by property A. Unfortunately, SOR techniques were not reliable in many of these problems, either converging very slowly or even diverging. Therefore, it is not a surprise that direct methods soon became the reference techniques in this field. The irregular sparsity pattern of the FE-discretized structural problems gave a renewed impulse to the development of more appropriate ordering strategies based on the graph theory and more efficient direct algorithms, leading in the 70s and early 80s to the formulation of the modern multifrontal solvers [13, 24, 25]. In the 80s the first pioneering commercial codes for FE structural computations became available and the solvers of choice were obviously direct methods because of their reliability and robustness. At this time the solution of sparse linear systems by direct techniques appears to be quite a mature field of research, well summarized in famous reference textbooks such as the ones by Duff et al. [26] and Davis [27].

In contrast, during the 60s and 70s iterative solvers were living their infancy. CG was rehabilitated as an iterative method by Reid in 1971 [28] who showed that for reasonably well-conditioned sparse SPD matrices the convergence could have been reached after far fewer than 𝑛 steps, being 𝑛 the size of 𝐴. This work set the pace for the extension of CG to nonsymmetric and indefinite problems, with the introduction of the earliest Krylov subspace methods such as the Minimal Residual (MINRES) [29] and the Biconjugate Gradient (Bi-CG) [30] algorithms. However, the crucial event for the future success of Krylov subspace methods was the publication in 1977 by Meijerink and van der Vorst of the Incomplete Cholesky CG (ICCG) algorithm [31]. Incomplete factorizations had been already introduced in the Soviet Union [32] and independently by Varga [33] in the late 50s, and the seminal ideas for improving the conditioning of the system matrix in the CG iteration can be found also in the original work by Hestenes and Stiefel [19] and in Engeli et al. [21], but Meijerink and van der Vorst were the first who put the things together and showed the great potential of the preconditioned CG. The idea was later extended by Kershaw [34] to an SPD matrix not necessarily of 𝑀-type and gained soon a great popularity. The key point was that CG as well as any other Krylov subspace method with a proper preconditioning could become competitive with the latest direct methods, though requiring much less memory and so being attractive for large three-dimensional (3D) simulations.

The 80s and 90s are the years of the great development of Krylov subspace methods. In 1986 Saad and Schultz introduced the Generalized Minimal Residual (GMRES) method [35] which soon became the algorithm of choice among the iterative methods for nonsymmetric linear systems. Some of the drawbacks of Bi-CG were addressed in 1989 by Sonnenveld who developed the Conjugate Gradient Squared (CGS) method [36], on the basis of which in 1992 van der Vorst presented the Biconjugate Gradient Stabilized (Bi-CGSTAB) algorithm [37]. Along with GMRES, Bi-CGSTAB is the most successful technique for nonsymmetric linear systems, with the advantage of relying on a short-term recurrence. As demonstrated in the famous Faber-Manteuffel theorem [38], Bi-CGSTAB is not optimal in the sense that it does not theoretically lead to an approximate solution with some minimum property in the current Krylov subspace; however, its convergence can be faster than GMRES. Another family of Krylov subspace methods includes the Quasi-Minimal Residual (QMR) algorithm [39], with its transpose-free (TFQMR) and symmetric (SQMR) variants [40, 41]. Several other variants of the methods above, including truncated, restarted, and flexible versions, along with nested optimal techniques, for example, [42, 43], have been introduced by several researchers until the late 90s. A recent review of Krylov subspace methods for linear systems is available in [44], while several textbooks have been written on this subject at the beginning of the 21st century, among which the most famous is perhaps the one by Saad [45].

Initially, Krylov subspace methods were seen with some suspicion by the practitioners, especially those coming from structural mechanics, because of their apparent lack of reliability. Numerical experiments in different fields, such as FE elastostatics, geomechanics, consolidation of porous media, fluid flow, and transport, for example, [46–51], showed a growing interest for these methods, mainly because they potentially could allow for the solution of much bigger and more detailed problems. Whereas direct methods typically scale poorly with the matrix size, especially in 3D models, Krylov subspace methods appeared to be virtually mandatory with large problems. On the other hand, it became soon clear that the key factor to improve the robustness and the computational efficiency of any iterative solver is preconditioning. Nowadays, it is quite well accepted that it is preconditioning, rather than the selected Krylov algorithm, the most important issue to address.

This is why research on the construction of effective preconditioners has significantly grown over the last two decades, while advances on Krylov subspace methods have progressively faded. Currently, preconditioning appears to be a much more active and promising research field than either direct and iterative solution methods, and particularly so within the context of the fast evolution of the hardware technology. On one hand, this is due to the understanding that there are virtually no limits to the available options for obtaining a good preconditioner. On the other hand, it is also generally recognized that an optimal general-purpose preconditioner is unlikely to exist, so that there is the possibility to improve the solver efficiency in different ways for any specific problem at hand within any specific computing environment. Generally, the knowledge of the governing physical processes, the structure of the resulting system matrix, and the available computer technology are factors that cannot be ignored in the design of an appropriate preconditioner. It is also recognized that theoretical results are few, and frequently “empirical” algorithms may work surprisingly well despite the lack of a rigorous foundation. This is why finding a good preconditioner for solving a sparse linear system can be viewed rather as “a combination of art and science” [45] than a rigorous mathematical exercise.

2. Basic Concepts

Roughly speaking, preconditioning means transforming system (1.1) into an equivalent mathematical problem which is expected to converge faster using an iterative solver. For instance, given a nonsingular matrix 𝑀, premultiplying (1.1) by 𝑀⁻¹ yields 𝑀⁻¹𝐴𝐱=𝑀⁻¹𝐛. (2.1) If the matrix 𝑀⁻¹𝐴 is better conditioned than 𝐴 for a Krylov subspace method, then 𝑀⁻¹ is the preconditioner and (2.1) denotes the left preconditioned system. Similarly, 𝑀⁻¹ can be applied on the right: 𝐴𝑀⁻¹𝐲=𝐛,𝐱=𝑀⁻¹𝐲, (2.2) thus producing a right preconditioned system. The use of left or right preconditioning depends on a number of factors and can produce quite different behaviors; see, for example, [45]. For instance, right preconditioning has the advantage that the residual of the preconditioned system is the same as the native system, while with left preconditioning it is not, and this can be important in the convergence check or using a residual minimization algorithm such as GMRES. By the way, with the right preconditioning a back transformation of the auxiliary variable 𝐲 must be carried out to resume the original unknown 𝐱. If the preconditioner can be written in a factorized form: 𝑀⁻¹=𝑀₂⁻¹𝑀₁⁻¹, (2.3) a split preconditioning can be also used: 𝑀₁⁻¹𝐴𝑀₂⁻¹𝐲=𝑀₁⁻¹𝐛,𝐱=𝑀₂⁻¹𝐲, (2.4) thus potentially exploiting the advantages of both left and right formulations.

Writing the preconditioned Krylov subspace algorithms is relatively straightforward. Simply, the basic algorithms can be implemented replacing 𝐴 with either 𝑀⁻¹𝐴, or 𝐴𝑀⁻¹, or 𝑀₁⁻¹𝐴𝑀₂⁻¹, and then back substituting the original variable 𝐱 where the auxiliary vector 𝐲 is used. It can be easily observed that it is not necessary to build the preconditioned matrix, which could be much less sparse than 𝐴, but just to compute the product of 𝑀⁻¹, or 𝑀₁⁻¹ and 𝑀₂⁻¹, if the factorized form (2.3) is used, by a vector. This operation is called application of the preconditioner.

Generally speaking, there are three basic requirements for obtaining a good preconditioner: (i)the preconditioned matrix should have a clustered eigenspectrum away from 0, (ii)the preconditioner should be as cheap to compute as possible, (iii)its application to a vector should be cost-effective. The origin of condition (i) relies on the convergence properties of CG. It is quite intuitive that if 𝑀⁻¹ resembles in some sense 𝐴⁻¹ the preconditioned matrix should be close to the identity, thus making the preconditioned system somewhat “easier” to solve and accelerating the convergence. If 𝐴 and 𝑀⁻¹ are SPD, it has been proved that the iteration count 𝑛_iter of the CG algorithm to go below the tolerance 𝜀 depends on the spectral conditioning number 𝜉 of the preconditioned matrix 𝐺, for example [52]: 𝑛_iter≤12𝜉(𝐺)2ln𝜀+1, (2.5) so that clustering the eigenvalues of 𝐺 away from 0 is of paramount importance to accelerate convergence. If 𝐴 is not SPD, things can be much more complicated. For example, the performance of GMRES depends also on the conditioning of the matrix of the eigenvectors of 𝐺, and the fact that the preconditioned matrix has a clustered eigenspectrum does not guarantee a fast convergence [53]. Anyway, it is quite a common experience that a clustered eigenspectrum away from 0 often yields a rapid convergence also in nonsymmetric problems.

The importance of conditions (ii) and (iii) depends on the specific problem at hand and may be highly influenced by the computer architecture. From a practical point of view, they can also be more important than condition (i). For example, if a sequence of linear systems has to be solved with the same matrix 𝐴 or with slight changes only, as it often may occur when using a Newton method for a set of nonlinear equations, the same preconditioner can be used several times, with the cost for its computation easily amortized. In this case, it could be convenient spending more time to build an effective preconditioner. On the other hand, the computational cost for the preconditioner application basically depends on its density. As 𝑀⁻¹ has to be applied once or twice per iteration, according to the selected Krylov subspace method, it is of paramount importance that the increased cost of each iteration is counterbalanced by an adequate reduction of their number. Typically, pursuing condition (i) reduces the iteration count but is in conflict with the need for a cheap computation and application of the preconditioner, so that in the end preconditioning efficiently the linear system (1.1) is always the result of an appropriate tradeoff.

An easy way to build a preconditioner is based on the decomposition of 𝐴. Recalling the definition of stationary iteration, 𝐱_𝑘+1=𝐱_𝑘+𝐶𝐫_𝑘, (2.6) where 𝐫_𝑘=𝐛−𝐴𝐱_𝑘 is the current residual, the matrix 𝐶 can be used as a preconditioner. Splitting 𝐴 as 𝐷−𝐸−𝐹, where 𝐷, −𝐸, and −𝐹 are the diagonal, the lower, and the upper parts of 𝐴, respectively, then 𝑀_𝐽⁻¹=𝐷⁻¹,𝑀_𝑆⁻¹=(𝐷−𝐸)⁻¹,𝑀_𝜔⁻¹=𝜔(𝐷−𝜔𝐸)⁻¹,𝑀_Ω⁻¹=𝜔(2−𝜔)(𝐷−𝜔𝐹)⁻¹𝐷(𝐷−𝜔𝐸)⁻¹ (2.7)(2.8)(2.9)(2.10) are known as the Jacobi, Gauss-Seidel, SOR, and Symmetric SOR preconditioners, respectively, where 𝜔 is the overrelaxation factor. Replacing 𝐸 with 𝐹 in (2.8) and (2.9) gives rise to the backward Gauss-Seidel and SOR preconditioners, respectively. The preconditioners obtained from the splitting of 𝐴 do not require additional memory and have a null computation cost, while their application can be simply done by a forward or backward substitution. With 𝑀_𝐽⁻¹ a simple diagonal scaling is actually required.

Because of their simplicity, Jacobi, Gauss-Seidel, SOR, and Symmetric SOR preconditioners are still used in some applications. Assume, for instance, that matrix 𝐴 has two clusters of eigenvalues due to the heterogeneity of the underlying physical problem. This is a quite common experience in geotechnical soil/structure or contact simulations, for example, [49, 54–58]. In this case the matrix 𝐴 can be written as 𝐴=𝐾𝐵₁𝐵₂𝐻, (2.11) where the square diagonal blocks are responsible of the two clusters of eigenvalues and usually arise in a natural way from the original problem. Factorizing 𝐴 in (2.11) according to Sylvester's theorem of inertia: 𝐵𝐴=𝐼0₂𝐾⁻¹𝐼𝐾00𝑆𝐼𝐾⁻¹𝐵₁0𝐼, (2.12) where 𝑆=𝐻−𝐵₂𝐾⁻¹𝐵₁ (2.13) is the Schur complement, a Generalized Jacobi (GJ) preconditioner [59], based on the seminal idea in [60], can be obtained as a diagonal approximation of the central block diagonal factor in (2.12): 𝑀⁻¹=diag(𝐾)00𝜑diag(𝑆)⁻¹, (2.14) where 𝜑 is a user-specified parameter and 𝑆 is an approximation of 𝑆 in (2.13) with diag(𝐾) used instead of 𝐾. Quite recently Bergamaschi and Martínez [61] have indirectly proved that an optimal value for 𝜑 is the ratio between the largest eigenvalues of 𝐾 and 𝑆. Another very simple algorithm which has proved effective in similar applications is based on SOR and Symmetric SOR preconditioners, for example, [62–64].

Though the algorithms above cannot compete with more sophisticated tools in terms of global performance [58, 65], they are an example of how even somewhat naive ideas can work and be useful for practitioners. Generally speaking, it is always a good idea to develop an effective preconditioner starting from the specific problem at hand. Interesting strategies have been developed approximating directly the differential operator the matrix 𝐴 arises from, by using a “nearby” PDE or a less accurate problem discretization, for example, [66–69]. For instance, popular preconditioners typically used within the transport community are often physically based, for example, [70–72].

For the developers of computer codes, however, using physically based preconditioners is not always desirable. In fact, it is usually simpler introducing the linear solver as a black box which is expected to work independently of the specific problem at hand, and especially so with less specialized codes such as the commercial ones. This is why a great interest has arisen also in the development of purely algebraic preconditioners that claim to be as much “general purpose” as possible. Algebraic preconditioners are usually robust algorithms which require the knowledge of the matrix 𝐴 only, independently of the underlying physical problem. This paper will concentrate on such a class of algorithms. Subdividing algebraic preconditioners into categories is not straightforward, as especially in recent years there have been several contaminations from adjacent fields with the development of “hybrid” approaches. Traditionally, and as already done by Benzi [2], algebraic preconditioners can be roughly divided into two classes, that is, incomplete factorizations and approximate inverses. The main distinction is that with the former approach 𝑀 is actually computed and 𝑀⁻¹ is applied but never formed explicitly, whereas with the latter 𝑀⁻¹ is directly built and applied. Within each class, moreover, different strategies can be pursued according to a number of factors, such as whether 𝐴 is SPD or not, its sparsity structure and size, the magnitude of its coefficients, and ultimately the computer hardware available to solve the problem, so that a large number of variants have been proposed by the researchers. As anticipated before, the unavoidable contamination from adjacent fields, for example, the direct solution methods, has made the boundaries between different groups of algorithms often vague and certainly questionable.

In this work, the traditional distinction between incomplete factorizations and approximate inverses is followed, describing the most successful algorithms belonging to each class and the most significant variants developed to address particular occurrences and increase efficiency. Then, two additional sections will be devoted to the most recent results obtained in the fields that currently appear to be the more active frontier of research in the area of algebraic preconditioning, that is, multigrid techniques and parallel algorithms. Finally, a few words will be spent on a special class of problems characterized by indefinite saddle-point matrices, which arise quite frequently in the applications and have attracted much interest in recent years.

3. Incomplete Factorizations

Given a nonsingular matrix 𝐴 such that 𝐴=𝐿𝑈, (3.1) the basic idea of incomplete factorizations is to approximate the triangular factors 𝐿 and 𝑈 in a cost-effective way. If 𝐿 and 𝑈 are computed so as to satisfy (3.1), for example, by a Gaussian elimination procedure, typically fill-in takes place, with 𝐿 and 𝑈 much less sparse than 𝐴. This is what limits the use of direct methods in large 3D problems. The approximations 𝐿 and 𝑈 of 𝐿 and 𝑈, respectively, can be obtained by simply discarding a number of fill-in entries according to some rules. The resulting preconditioner 𝑀⁻¹=𝑈⁻¹𝐿⁻¹ (3.2) is generally denoted as Incomplete LU (ILU). Quite obviously, 𝑀⁻¹ in (3.2) is never built explicitly as it is much denser than both 𝐿 and 𝑈. Its application to a vector is therefore performed by forward and backward substitutions. Moreover, the ILU factorized form (3.2) is well suited for split preconditioning. In case 𝐴 is SPD, the upper incomplete factor 𝑈 is replaced by 𝐿^𝑇 and the related preconditioner is known as Incomplete Cholesky (IC).

The native ILU algorithm runs as follows. Define a set 𝒮 of position (𝑖,𝑗), with 1≤𝑖,𝑗≤𝑛, where either 𝐿 if 𝑖>𝑗 or 𝑈 if 𝑗>𝑖 has a nonzero entry. The set 𝒮 is also denoted as the nonzero pattern of the incomplete factors. To avoid breakdowns in the ILU computation, it is generally recommended to include all the diagonal entries in 𝒮. Then, make a copy of 𝐴 over an auxiliary matrix which will contain 𝐿 and 𝑈 in the lower and upper parts, respectively, at the end of the computation, and set to zero all entries 𝑎_𝑖𝑗 such that (𝑖,𝑗)∉𝒮. If the main diagonal belongs to 𝐿, then it is implicitly assumed that the diagonal of 𝑈 is unitary and vice versa. For every position (𝑖,𝑗)∈𝒮, the following update is computed: 𝑎_𝑖𝑗⟵𝑎_𝑖𝑗−𝑎_𝑖𝑘𝑎⁻¹_𝑘𝑘𝑎_𝑘𝑗 (3.3) with 𝑘=1,…,𝑖−1. At the end of the loop, the ILU factors are stored in the copy of 𝐴.

The procedure described above is clearly an incomplete Gauss elimination. Hence, no surprise that several implementation tricks developed to improve the efficiency of direct solution methods can be borrowed in order to reduce the computational cost of the ILU computation. For example, the algorithm can follow either the KIJ, or the IKJ, or the IJK elimination variants. To give an idea, Algorithm 1 provides a sketch of the sequence of operations required by the IKJ implementation, which is generally preferred whenever the sparse matrix 𝐴 along with both 𝐿 and 𝑈 is stored in a compact form by rows. It can be easily observed that Algorithm 1 proceeds by rows, computing first the row of 𝐿 and then that of 𝑈 and accessing the previous rows only to gather the required information (Figure 1). Many other efficient variants have been developed, for example, storing 𝐴 in a sparse skyline format or computing 𝐿 by columns [73, 74].

Input: Matrix 𝐴, the nonzero pattern 𝒮 Output: Matrix 𝐴 containing 𝐿 and 𝑈 for  each(𝑖,𝑗)∉𝒮  do   𝑎𝑖𝑗=0 end for  𝑖=2,...,𝑛  do   for  𝑘=1,…,𝑖−1 and (𝑖,𝑘)∈𝒮  do   𝑎𝑖𝑘←𝑎𝑖𝑘/𝑎𝑘𝑘   for  𝑗=𝑘+1,...,𝑛and(𝑖,𝑗)∈𝒮  do    𝑎𝑖𝑗←𝑎𝑖𝑗−𝑎𝑖𝑘𝑎𝑘𝑗   end  end end

Algorithm 1 

IKJ variant of ILU factorization with a static pattern 𝒮.

Figure 1 

IKJ variant of ILU factorization.

3.1. Fill-In Strategies

The several existing versions of the ILU preconditioner basically differ for the rules followed to select the retained entries in 𝐿 and 𝑈. The simplest idea is to define 𝒮 statically at the beginning of the algorithm and equal to the set of the nonzero positions in 𝐴. This idea was originally implemented by Meijerink and van der Vorst [31] in their 1977 seminal work, giving rise to what is traditionally referred to as ILU(0), or IC(0) for SPD matrices. This simple choice is easy to implement and very cost-effective, as the preconditioner construction is quite cheap and its application to a vector costs as much as a matrix-by-vector product with 𝐴. Moreover, ILU(0) can be quite efficient, especially for matrices arising from the discretization of elliptic problems, for example, [75–77].

Unfortunately, there are a number of problems where the ILU(0) preconditioner converges very slowly or even fails. There are several different reasons for this, depending on both the unstable computation and the unstable application of 𝐿 and 𝑈, as was experimentally found by Chow and Saad on a large number of test matrices [78]. For example, ILU(0) often exhibits a bad behavior with structural matrices and the highly nonsymmetric and indefinite matrices arising from computational fluid dynamics. In these cases the nonzero pattern of 𝐴 proves a too crude approximation of the exact 𝐿 and 𝑈 structure, so that a possible remedy is to allow for a larger number of entries in 𝐿 and 𝑈 in order to make them closer to 𝐿 and 𝑈. In the limiting case, the exact factors are computed and the solution to (1.1) is obtained in just one iteration. Obviously, enlarging the ILU pattern on one hand reduces the iteration count but on the other increases the cost of computation and application of the preconditioner, so an appropriate tradeoff must be found.

There are several different ways for enlarging the initial pattern 𝒮, or computing it in a dynamic way. In any specific problem, the winning algorithm is the one that proves able to capture the most significant entries in 𝐿 and 𝑈 at the lowest cost. One of the first attempts for enlarging dynamically 𝒮 is based on the so-called level-of-fill concept [79, 80]. The idea is to assign an integer value ℓ_𝑖𝑗, denoted as level-of-fill, to each entry in position (𝑖,𝑗) computed during the ILU construction process. The lower the level, the more important the entry. The initial level is set to ℓ_𝑖𝑗=0if𝑎_𝑖𝑗≠0∞if𝑎_𝑖𝑗=0. (3.4) During the ILU construction, with reference for instance to the Algorithm 1, the level-of-fill of each entry is updated according to the level-of-fill of the entries used for its computation using the following rule: ℓ_𝑖𝑗⟵minℓ_𝑖𝑗,ℓ_𝑖𝑘+ℓ_𝑘𝑗+1. (3.5) In practice, any position corresponding to a nonzero entry in 𝐴 remains with a zero level-of-fill. The terms computed using entries with a zero level-of-fill have a level equal to 1. Those built using one entry of the first level have level 2 and so on. The new level-of-fill of each entry is computed adding 1 to the sum of the levels of the entries used for it. This kind of algorithm creates a hierarchy among the potential entries that allows for selecting the “most important” nonzero terms only. The resulting preconditioner is denoted as ILU(𝑝), where 𝑝>0 is the integer user-specified level-of-fill of the entries retained into each row of 𝐿 and 𝑈. It follows immediately that ILU(0) coincides exactly with the zero fill-in ILU previously defined.

In most applications, ILU(1) is already a significant improvement to ILU(0). The experience shows that the fill-in increases quite rapidly with 𝑝, so that for 𝑝≥2 the cost for building and applying the preconditioner can become a price too high to pay. Moreover it can be observed that at each level there are also many terms which are very small in magnitude, thus providing a small contribution to the preconditioning of the original matrix 𝐴. Hence, a natural way to avoid such an occurrence is to add a dropping rule according to which the terms below a user-specified tolerance 𝜏 are neglected. Several criteria have been defined for selecting 𝜏 and the way the ILU entries are discarded. A usual choice is to compute the 𝑖th row of both 𝐿 and 𝑈 and drop an entry whenever smaller than 𝜏‖𝐚_𝑖‖, being ‖⋅‖ an appropriate vector norm and 𝐚_𝑖 the 𝑖th row of 𝐴, but other strategies have been also attempted, such as comparing the entry 𝑙_𝑖𝑗 of 𝐿 with 𝜏𝑙_𝑖𝑖𝑙_𝑗𝑗 [81] or with an estimate of the norm of the 𝑗th column of 𝐿⁻¹ [82, 83].

Similarly to the level-of-fill approach, a drawback of the drop tolerance strategy is that the amount of fill-in of 𝐿 and 𝑈 is strongly dependent on the selected user-specified parameter and the specific problem at hand and is unpredictable a priori. To address such an issue a second user-specified parameter 𝜌 can be introduced, denoting the largest number of entries that can be retained per row. In some versions, 𝜌 is also defined as the number of entries added at most per row in excess to the number of entries of the same row of 𝐴. The use of this dual threshold strategy, which may be quite simple and effective, has been proposed by Saad [84], and the related preconditioner is referred to as ILUT(𝜌,𝜏). A sketch of the sequence of steps required for its computation is provided in Algorithm 2.

Input: Matrix 𝐴, the user-specified parameters 𝜌 and 𝜏 Output: Matrices 𝐿 and 𝑈 for  𝑖=1,…,𝑛  do  𝐰=𝐚𝑖  for  𝑘=1,…,𝑖−1 and 𝑤𝑘≠0  do   𝑤𝑘←𝑤𝑘/𝑎𝑘𝑘   if  𝑤𝑘≥𝜏‖𝐚𝑖‖2  then     𝐰←𝐰−𝑤𝑘𝐮𝑘   else    𝑤𝑘=0   end  end  Retain the 𝜌 largest entries in  w  for  𝑗=1,…,𝑖−1  do    𝑙𝑖𝑗=𝑤𝑗  end  for  𝑗=𝑖,…,𝑛  do     𝑢𝑖𝑗=𝑤𝑗  end end

Algorithm 2 

ILUT(𝜌,𝜏) computation.

As the drop tolerance may be quite sensitive to the specific problem and so difficult to implement in a black box solver, ILU variants that use only the fill-in parameter 𝜌 have been proposed [85, 86]. These versions are generally outperformed by ILUT(𝜌,𝜏) but are much easier to set and generally less sensitive to the user-specified parameter, allowing the user for a thorough control of the memory occupation. More recent ILU variants developed with the aim of improving the preconditioner performance and better exploiting the hardware potential are based on multilevel techniques, dense block identification, and adaptive selections of the setup parameters [87–89], with significant contaminations from the methods developed for the latest multifrontal solvers.

3.2. Stabilization Techniques

Both Algorithms 1 and 2 require the execution of divisions where the denominator, that is, the pivot, is not generally guaranteed to be nonzero. It is well known that in finite arithmetics also small pivots can create numerical difficulties. If 𝐴 is SPD, the pivot has to be also positive in order to produce a real incomplete factorization. These issues are very important because they can jeopardize the existence of the preconditioner itself and so the robustness of the overall iterative algorithm.

Meijerink and van der Vorst [31] demonstrated that for an M-matrix a zero fill-in IC is guaranteed to exist. However, if 𝐴 is SPD but not of M-type zero or negative pivots may arise, thus theoretically preventing from the IC computation. Typically, such an inconvenience may be avoided by increasing properly the fill-in. In this way the incomplete factor tends to the exact factor, which is guaranteed to be SPD. However, this may generate a fill-in degree that is unacceptably high, lowering the performance of the preconditioned solver. To avoid this undesirable occurrence a number of tricks have been proposed which can greatly increase the IC robustness. The first, and maybe simplest, idea was advanced by Kershaw [34] who just suggested to replace zero or negative pivots with an arbitrary positive value, for example, the last valid pivot. The implementation of this trick is straightforward, but it seldom provides good results as the resulting IC often exhibits a very low quality. In fact, recall from Algorithm 1 that modifying arbitrarily a pivot in the 𝑘th row means influencing the computation of all the following rows. Another simple strategy recently advanced in [90] relies on avoiding the computation of the square root of the pivot, required in the IC algorithm, thus eliminating a source of breakdown. In [91] a diagonal compensated reduction of the positive off-diagonal entries is suggested in order to transform the native matrix into an M-matrix before the incomplete factorization process is performed. A different algorithm is proposed in [92] with the incomplete factors obtained through an 𝐴-orthogonalization instead of a Cholesky elimination.

A famous stabilization technique to ensure the IC existence is the one advanced by Ajiz and Jennings [93], that has proved quite effective especially in structural applications [50, 90]. Following this approach, we can write 𝐴=𝐿𝐿^𝑇=𝐿𝐿^𝑇+𝑅, (3.6) where 𝑅 is a residual matrix collecting with all the entries discarded during the incomplete procedure. From (3.6), 𝐿𝐿^𝑇 is actually the exact factorization of 𝐴−𝑅. As 𝐴 is SPD, it follows immediately that a sufficient condition for 𝐿 to be real is that 𝑅 is negative semidefinite. Suppose that the first column of 𝐿 has been computed and the term 𝑙_𝑗1=𝑎_𝑗1/𝑎₁₁ has to be dropped along with the corresponding 𝑙_1𝑗 term in 𝐿^𝑇. Hence, the incomplete factorization is equivalent to the exact factorization of 𝐴 where the entries 𝑎_1𝑗 and 𝑎_𝑗1 are replaced by zero, and the error matrix 𝑅 contains 𝑎_𝑗1 in position 1𝑗 and 𝑗1. For example, in a 3×3 matrix with 𝑗=2, we have 𝑎𝐴=¹¹𝑎₁₂𝑎₁₃𝑎₂₁𝑎₂₂𝑎₂₃𝑎₃₁𝑎₃₂𝑎₃₃=𝑎₁₁0𝑎₁₃0𝑎₂₂𝑎₂₃𝑎₃₁𝑎₃₂𝑎₃₃+0𝑎₁₂0𝑎₂₁00=000𝐴+𝑅. (3.7) The residual matrix 𝑅 can be forced to be negative semidefinite by adding two coefficients 𝛼 and 𝛽 in position (1,1) and (𝑗,𝑗), such that the submatrix 𝛼𝑎_1𝑗𝑎_𝑗1𝛽 (3.8) is negative semidefinite. The matrix 𝐴 must be modified accordingly by subtracting 𝛼 and 𝛽 from the corresponding diagonal terms. Though different choices are possible, for example, [50, 93], the option 𝛼=𝛽=−|𝑎_1𝑗| is particularly attractive because in this way the sum of the absolute values of the arbitrarily introduced entries |𝛼|+|𝛽| is minimized. This procedure ensures that 𝐿𝐿^𝑇 is the exact factorization of a positive definite matrix, hence its existence is guaranteed, but obviously such a matrix can be quite different from 𝐴. The quality of this stabilization strategy basically depends on the quality of 𝐴 as an approximation of 𝐴.

Unstable factors and negative pivots typically arise if the diagonal terms of 𝐴 are remarkably smaller than the off-diagonal ones, or if there are significant contrasts in the coefficient magnitude. The latter occurrence is particularly frequent in heterogeneous problems or in models coupling different physical processes, such as multiphase flow in deformable media, coupled flow and transport, themomechanical models, fluid-structure interactions, or multibody systems, for example, [94–97]. All these problems give naturally rise to a multilevel matrix: 𝐾𝐴=¹𝐵₁₂𝐵₁₃⋯𝐵_1ℓ𝐵₂₁𝐾₂𝐵₂₃⋯𝐵_2ℓ𝐵₃₁𝐵₃₂𝐾₃⋯𝐵_3ℓ𝐵⋮⋮⋮⋱⋮_ℓ1𝐵_ℓ2𝐵_ℓ3⋯𝐾_ℓ, (3.9) where the system unknowns can be subdivided into contiguous groups. By distinction with multigrid methods, a multilevel approach assumes that the different levels do not necessarily correspond to model subgrids. Several strategies have been proposed to handle properly the matrix levels, for example, [98–102]. Generally speaking, all such algorithms can be viewed as approximate Schur complement methods that differ for the strategy used to perform the level subdivision and the restriction and prolongation operators adopted [103, 104]. For a recent review, see, for instance, [105]. The basic idea of multilevel ILU proceeds as follows. Consider the partial factorization of 𝐴 in (3.9), where 𝐴₁ is the submatrix obtained from 𝐴 without the first level, and 𝐵₁₂ and 𝐵₂₁ are the rectangular submatrices collecting the blocks 𝐵_1𝑗 and 𝐵_𝑗1, 𝑗=2,…,ℓ: 𝐾𝐴=¹𝐵₁₂𝐵₂₁𝐴₁=𝐿₁0𝐵₂₁𝑈₁⁻¹𝐼𝐼00𝑆₁𝑈₁𝐿₁⁻¹𝐵₁₂0𝐼. (3.10) In (3.10) 𝐾₁=𝐿₁𝑈₁ and 𝑆₁ is the Schur complement: 𝑆₁=𝐴₁−𝐵₂₁𝑈₁⁻¹𝐿₁⁻¹𝐵₁₂. (3.11) Replacing 𝐿₁ and 𝑈₁ with the ILU factors 𝐿₁ and 𝑈₁ of 𝐾₁ gives rise to a partial incomplete factorization of 𝐴 that can be used as a preconditioner: 𝑀⁻¹=𝑈₁𝐻₁₂0𝐼⁻¹0𝑆𝐼0₁₋₁𝐿₁0𝐻₂₁𝐼⁻¹, (3.12) where 𝐻₁₂=𝐿₁⁻¹𝐵₁₂𝐻₂₁=𝐵₂₁𝑈₁⁻¹,𝑆₁=𝐴₁−𝐻₂₁𝐻₁₂. (3.13)(3.14) As 𝐻₁₂ and 𝐻₂₁ tend to be much denser than 𝐵₁₂ and 𝐵₂₁, some degree of dropping can be conveniently enforced, thus replacing these blocks with 𝐻₁₂ and 𝐻₂₁, respectively. Recalling the level structure of 𝐴₁, the inverse of the approximated Schur complement 𝑆₁ can be also replaced by a partial incomplete factorization with the same structure as 𝑀⁻¹ in (3.12). Repeating recursively the procedure above starting from 𝑆₀=𝐴, at the 𝑖th level, 𝑆_𝑖⁻¹ is replaced by 𝑆_𝑖⁻¹≃𝑈_𝑖+1𝐻_{𝑖+1,𝑖+2}0𝐼⁻¹0𝑆𝐼0_−1𝑖+1𝐿_𝑖+10𝐻_{𝑖+2,𝑖+1}𝐼⁻¹, (3.15) thus giving rise to an (𝑖+1)th level Schur complement 𝑆_𝑖+1=𝐴_𝑖+1−𝐻_{𝑖+2,𝑖+1}𝐻_{𝑖+1,𝑖+2}. The algorithm stops when the size of 𝑆_𝑖+1 equals 0, that is, when 𝑖+1 equals the user specified number of levels ℓ. Multilevel ILU preconditioners are generally much more robust than standard ILU and can be very effective provided that a suitable level subdivision is defined. If 𝐴 is SPD, however, it is necessary to enforce the positive definiteness of all approximate Schur complements, which is no longer ensured when incomplete decompositions and dropping procedures are used. Recently, Janna et al. [106] have developed a stabilization strategy which can guarantee the positive definiteness of each 𝑆_𝑖. Basically, they extend to levels of the factorization procedure introduced by Tismenetsky [107] where the Schur complement is updated by a quadratic formula. The theoretical robustness of Tismenetsky's algorithm has been proved in [108]. Recalling (3.6), the matrix 𝑆_𝑖 can be additively decomposed as 𝑆_𝑖=𝑆_𝑖,11𝑆_𝑖,12𝑆^𝑇_𝑖,12𝑆_𝑖,22=𝑆_𝑖+𝑅_𝑖=𝑆_𝑖,11𝑆_𝑖,12𝑆^𝑇_𝑖,12𝑆_𝑖,22+𝑅000, (3.16) where 𝑅 is the residual matrix obtained from the IC decomposition of 𝑆_𝑖,11. Assuming that an appropriate stabilization technique has been used, for example, the one by Ajiz and Jennings [93], 𝑅_𝑖 is negative semidefinite, so 𝑆_𝑖 is positive definite. Let us collect all the entries dropped in the 𝐻_{𝑖+1,𝑖+2} computation in the error matrix 𝐸_𝐻, thus writing 𝐻_{𝑖+1,𝑖+2}=𝐻_{𝑖+1,𝑖+2}+𝐸_𝐻. The Schur complement of 𝑆_𝑖 can be now rewritten as 𝑆_𝑖+1=𝑆_𝑖,22−𝐻^𝑇_{𝑖+1,𝑖+2}𝐻_{𝑖+1,𝑖+2}−𝐸^𝑇_𝐻𝐻_{𝑖+1,𝑖+2}−𝐻^𝑇_{𝑖+1,𝑖+2}𝐸_𝐻−𝐸^𝑇_𝐻𝐸_𝐻. (3.17) Ignoring the last three terms of (3.17) we obtain 𝑆_𝑖+1=𝑆_𝑖,22−𝐻^𝑇_{𝑖+1,𝑖+2}𝐻_{𝑖+1,𝑖+2}, (3.18) that is the standard Schur complement usually computed in a multilevel algorithm; see (3.14). As 𝐸^𝑇_𝐻𝐻_{𝑖+1,𝑖+2} is generally indefinite, 𝑆_𝑖+1 in (3.18) may be indefinite too. As a remedy, add and subtract 𝐻^𝑇_{𝑖+1,𝑖+2}𝐻_{𝑖+1,𝑖+2} to the right-hand side of (3.17): 𝑆_𝑖+1=𝑆_𝑖,22+𝐻^𝑇_{𝑖+1,𝑖+2}𝐻_{𝑖+1,𝑖+2}−𝐻_{𝑖+1,𝑖+2}+𝐸_𝐻^𝑇𝐻_{𝑖+1,𝑖+2}−𝐻^𝑇_{𝑖+1,𝑖+2}𝐻_{𝑖+1,𝑖+2}+𝐸_𝐻−𝐸^𝑇_𝐻𝐸_𝐻. (3.19) Recalling that 𝐻_{𝑖+1,𝑖+2}+𝐸_𝐻=𝐻_{𝑖+1,𝑖+2}=𝐿⁻¹_𝑖+1𝑆_𝑖,12, (3.20) (3.19) becomes 𝑆_𝑖+1=𝑆_𝑖,22+𝐻^𝑇_{𝑖+1,𝑖+2}𝐻_{𝑖+1,𝑖+2}−𝑆^𝑇_𝑖,12𝐿^−𝑇_𝑖+1𝐻_{𝑖+1,𝑖+2}−𝐻^𝑇_{𝑖+1,𝑖+2}𝐿⁻¹_𝑖+1𝑆_𝑖,12−𝐸^𝑇_𝐻𝐸_𝐻. (3.21) Neglecting the last term in (3.21) we obtain another expression for the Schur complement: 𝑆_𝑖+1=𝑆_𝑖,22+𝐻^𝑇_{𝑖+1,𝑖+2}𝐻_{𝑖+1,𝑖+2}−𝑆^𝑇_𝑖,12𝐿^−𝑇_𝑖+1𝐻_{𝑖+1,𝑖+2}−𝐻^𝑇_{𝑖+1,𝑖+2}𝐿⁻¹_𝑖+1𝑆_𝑖,12, (3.22) which is always SPD. In fact, (3.21) yields 𝑆_𝑖+1=𝑆_𝑖+1+𝐸^𝑇_𝐻𝐸_𝐻, (3.23) which is SPD independently of the degree of dropping enforced on 𝐻_{𝑖+1,𝑖+2}. The procedure above is generally more expensive than the standard one but is also more robust.

If 𝐴 is not SPD, enforcing the positivity of pivots or the positive definiteness of the Schur complements is not necessary. However, the experience shows that, in many problems, especially those with a high degree of nonsymmetry and a lack of diagonal dominance, ILU can be much more unstable in both the computation and the application, leading to frequent solver failures [78]. In this case no rigorous methods have been introduced to avoid numerical instabilities. The effect of a preliminary reordering and/or scaling of 𝐴 has been debated for quite a long time, without achieving a shared conclusion. The fact that the matrix nonzero pattern has an effect on the quality of the resulting ILU preconditioner, similarly to what happens with direct methods, is quite well understood, but not all the methods performing well with direct techniques appear to be efficient for ILU, for example, [90, 109–112]. A similar result is obtained with a preliminary scaling. As proved in [113], scaling the native matrix does not change the spectral properties of an ILU preconditioned matrix; however it can greatly stabilize the ILU computation by reducing significantly the impact of the round-off errors. Useful scaling strategies can be taken from the GJ preconditioning of (2.14) or the Least Square Log algorithm introduced in [18].

4. Approximate Inverses

Incomplete factorizations can be very powerful preconditioners; however the issues concerning their existence and numerical stability can undermine their robustness in several examples. One reason for such a weakness can be understood using the following simple argument. Consider the incomplete factorization 𝐿𝑈 along with the corresponding residual matrix for a general matrix 𝐴: 𝐴=𝐿𝑈+𝑅. (4.1) Left and right multiplying both sides of (4.1) by 𝐿⁻¹ and 𝑈⁻¹, respectively, yields 𝐿⁻¹𝐴𝑈⁻¹=𝐼+𝐿⁻¹𝑅𝑈⁻¹. (4.2) The matrix controlling the convergence of a preconditioned Krylov subspace method is actually 𝐿⁻¹𝐴𝑈⁻¹ which differs from the identity by 𝐿⁻¹𝑅𝑈⁻¹. This means that a residual matrix close to 0 is not sufficient to guarantee a fast convergence. In fact, if the norms of either 𝐿⁻¹ or 𝑈⁻¹ are large, 𝐿⁻¹𝐴𝑈⁻¹ can be anyway far from the identity and yields a slow convergence, that is, even an accurate incomplete factorization can provide a poor outcome. This is why problems with strongly nonsymmetric matrices and a lack of diagonal dominance incomplete factorizations often have a bad reputation.

To cope with these drawbacks the researchers in the past have developed a second big class of preconditioners with the aim of computing an explicit form for 𝑀⁻¹, thus avoiding the solution of the linear system needed to apply the preconditioner. Generally speaking, an approximate inverse is an explicit approximation of 𝐴⁻¹. Quite obviously, it has to be sparse in order to maintain a workable cost for its computation and application. It is no surprise that several different methods have been proposed to build an approximate inverse. The most successful ones can be roughly classified according to whether 𝑀⁻¹ is computed monolithically or as the product of two matrices. In the first case, the approximate inverse is generally computed by minimizing the Frobenius norm: 𝐼−𝐴𝑀⁻¹_𝐹or𝐼−𝑀⁻¹𝐴_𝐹, (4.3) thus obtaining either a right or a left approximate inverse, which is generally different. In the second case, the basic idea is to find an explicit approximation of the inverse of the triangular factors of 𝐴: 𝑀₁⁻¹≃𝐿⁻¹,𝑀₂⁻¹≃𝑈⁻¹, (4.4) and then use 𝑀⁻¹ in the factored form (2.3). Different approaches have been advanced within each class and also in between these two groups.

Early algorithms for computing an approximate inverse via the Frobenius norm minimization appeared as far back as in the 70s [114, 115]. The basic idea is to select a nonzero pattern 𝒮 for 𝑀⁻¹ including all the positions where nonzero terms may be expected. Denoting by 𝒲_𝑆 the set of matrices with nonzero pattern 𝒮, any algorithm belonging to this class looks for a matrix 𝑀⁻¹∈𝒲_𝑆 such that ‖𝐼−𝐴𝑀⁻¹‖_𝐹 is minimum. Recalling the definition of Frobenius norm: ‖𝐴‖_𝐹=^𝑛𝑛_𝑗=1𝑖=1𝑎²_𝑖𝑗, (4.5) it follows that 𝐼−𝐴𝑀⁻¹²_𝐹=^𝑛_𝑗=1𝐞_𝑗−𝐴𝐦_𝑗²₂, (4.6) where 𝐞_𝑗 and 𝐦_𝑗 are the 𝑗th column of 𝐼 and 𝑀⁻¹, respectively. The sum at the right-hand side of (4.6) is minimum only if each addendum is minimum. Therefore, the computation of 𝑀⁻¹ is equivalent to the solution of 𝑛 least-squares problems with the constraint that 𝐦_𝑗 has nonzero components only in the positions 𝑖 such that (𝑖,𝑗)∈𝒮. This requirement reduces enormously the computational cost required to build 𝑀⁻¹. Denote by 𝒥 the set 𝒥=𝑖∶𝑖,𝑗∈𝒮 (4.7) and 𝐦_𝑗[𝒥] the subvector of 𝐦_𝑗 made of the components included in 𝒥. In the product 𝐴𝐦_𝑗 only the columns of 𝐴 with index in 𝒥 will provide a nonzero contribution. Moreover, as 𝐴 is sparse, only a few rows will have a nonzero term in at least one of these columns (Figure 2). Denoting by ℐ the set ℐ=𝑖∶𝑎_𝑖𝑗≠0with𝑗∈𝒥 (4.8) and 𝐴[ℐ,𝒥] the submatrix of 𝐴 made of the coefficients 𝑎_𝑖𝑗 such that 𝑖∈ℐ and 𝑗∈𝒥, each least square problem at the right-hand side of (4.6) reads 𝑒_𝑗[ℐ]−𝐴[ℐ,𝒥]𝐦_𝑗[𝒥]²₂→min (4.9) and is typically quite small. The procedure illustrated above allows for the computation of a right approximate inverse. To build a left approximate inverse just observe that 𝐼−𝑀⁻¹𝐴_𝐹=𝐼−𝐴^𝑇𝑀^−𝑇_𝐹. (4.10) Hence, (4.9) can still be used for any 𝑗 replacing 𝐴 with 𝐴^𝑇 and obtaining the rows of 𝑀⁻¹. It turns out that for strongly nonsymmetric matrices the right and left approximate inverse can be quite different.

Figure 2 

Schematic representation of a matrix-vector product subject to sparsity constraints.

The main difference in the algorithms developed within this class relies on how the nonzero pattern 𝒮 is selected, and this is indeed the key factor for the overall quality of the approximate inverse. A more detailed discussion on this point will be provided in the next section. Some of the most popular algorithms in this class have been advanced in [116–121], each one with its pros and cons. Probably the most successful approach is the one proposed by Grote and Huckle [117] which is known as SPAI (SParse Approximate Inverse). Basically, the SPAI algorithm follows the lines described above dynamically enlarging 𝒥 for each column of 𝑀⁻¹ until the least square problem (4.9) is solved with a prescribed accuracy. SPAI has proved to be quite a robust method allowing for a good control of its quality by the user. Unfortunately, its computation can be quite expensive also on a parallel machine. A sketch of the SPAI construction is reported in Algorithm 3.

Input: Matrix 𝐴, the initial nonzero pattern 𝒮 and the tolerance 𝜀 Output: Matrix 𝑀−1 for  𝑗=1,…,𝑛  do  Compute 𝒥={𝑖∶(𝑖,𝑗)∈𝒮}  Compute ℐ={𝑖∶𝑎𝑖𝑘≠0,𝑘∈𝒥}  Gather 𝐴[ℐ,𝒥] from 𝐴  Solve the least square problem ‖𝐞𝑗[ℐ]−𝐴[ℐ,𝒥]𝐦𝑗[𝒥]‖2→min  𝐫𝑗=𝐞𝑗[ℐ]−𝐴[ℐ,𝒥]𝐦𝑗[𝒥]  while  ‖𝐫𝑗‖2>𝜀  do   Enlarge 𝒥 and update ℐ   Gather 𝐴[ℐ,𝒥] from 𝐴   Solve the least square problem ‖𝐞𝑗[ℐ]−𝐴[ℐ,𝒥]𝐦𝑗[𝒥]‖2→min   𝐫𝑗=𝐞𝑗[ℐ]−𝐴[ℐ,𝒥]𝐦𝑗[𝒥]  end end

Algorithm 3 

Computation of the SPAI preconditioner.

A different way to build an approximate inverse of 𝐴 relies on approximating the inverse of its triangular factors. Assuming that all the leading principal minors of 𝐴 are not null, that is, 𝐴 can be factorized as the product 𝐿𝐷𝑈 where 𝐿 is unit lower triangular, 𝐷 is diagonal, and 𝑈 is upper unit triangular, then 𝐿⁻¹𝐴𝑈⁻¹=𝐷 (4.11) implying that the rows of 𝐿⁻¹ and the columns of 𝑈⁻¹ are two sets of 𝐴-biconjugate vectors. In other words, denoting by 𝐰_𝑖 and 𝐳_𝑗 the 𝑖th row of 𝐿⁻¹ and the 𝑗th column of 𝑈⁻¹, respectively, (4.11) means that 𝐰^𝑇_𝑖𝐴𝐳_𝑗=0∀𝑖,𝑗=1,…,𝑛,𝑖≠𝑗. (4.12) If 𝑊 and 𝑍 are the matrices collecting all vectors 𝐰_𝑖 and 𝐳_𝑗, 𝑖,𝑗=1,…,𝑛, satisfying the condition (4.12) and 𝐷=𝑊^𝑇𝐴𝑍, then the inverse of 𝐴 reads 𝐴⁻¹=𝑍𝐷⁻¹𝑊^𝑇=^𝑛_𝑘=1𝐳_𝑘𝐰^𝑇_𝑘𝑑_𝑘 (4.13) and can be computed explicitly. There are infinite sets of vectors satisfying (4.12), and they can be computed by means of a biconjugation process, for example, via a two-sided Gram-Schmidt algorithm, starting from two arbitrary initial sets. Quite obviously, these vectors are dense and their explicit computation is practically unfeasible. However, if the process is carried out incompletely, for instance, dropping the entries of 𝐰_𝑘 and 𝐳_𝑘, 𝑘=1,…,𝑛, whenever smaller than a prescribed tolerance, the matrices 𝑊 and 𝑍 can be acceptable approximations of 𝐿^−𝑇 and 𝑈⁻¹ though preserving a workable sparsity. The resulting preconditioner is denoted as AINV (Approximate INVerse) and has been introduced in the late 90s by Benzi and coauthors [122–124]. In the following years the native algorithm has been further developed by other researchers as well, for example, [125–130]. A sketch of the basic procedure to build the AINV preconditioner is provided in Algorithm 4.

Input: Matrix 𝐴, the tolerance 𝜏 Output: Matrices 𝑊, 𝑍 and 𝐷 such that 𝑀−1=𝑍𝐷−1𝑊𝑇 Set p  =  q  =  0, 𝑍=𝑊=𝐼 for   𝑘=1,…,𝑛  do  𝑝𝑘=𝐰𝑇𝑘𝐴𝐳𝑘,𝑞𝑘=𝐰𝑇𝑘𝐴𝐳𝑘  for  𝑖=𝑘+1,…,𝑛  do   𝑝𝑖=(𝐰𝑇𝑖𝐴𝐳𝑘)/𝑝𝑘,𝑞𝑖=(𝐰𝑇𝑘𝐴𝐳𝑖)/𝑞𝑘   𝐰𝑖←𝐰𝑖−𝑝𝑖𝐰𝑘,𝐳𝑖←𝐳𝑖−𝑞𝑖𝐳𝑘   for  𝑗=1,…,𝑖  do    Drop the 𝑗th component of 𝐰𝑖 and 𝐳𝑖 if smaller than 𝜏    end   end end for  𝑘=1,…,𝑛  do    𝑑𝑘=𝑝𝑘 end

Algorithm 4 

Computation of the AINV preconditioner.

Quite obviously, if 𝐴 is SPD, then 𝑊=𝑍, and in Algorithm 4 only the update of 𝐳_𝑖 must be carried out. The resulting AINV preconditioner can be written as 𝑀⁻¹=𝑍𝑍^𝑇, is SPD, and exists if 𝑝_𝑘≠0 for all 𝑘=1,…,𝑛. Unfortunately, the same does not hold true for nonfactored sparse approximate inverses based on minimizing the Frobenius norm in (4.3). In general, using an algorithm developed along the lines of SPAI there is no guarantee that 𝑀⁻¹ is SPD if 𝐴 is so. This is quite a strong limitation, as it prevents from using the powerful CG method to solve an SPD linear system. Such a weak point motivated the development of the Factorized Sparse Approximate Inverse (FSAI) preconditioner by Kolotilina and Yeremin [131]. The FSAI algorithm lies somewhat between the two classes of approximate inverses previously identified, as it produces a factored preconditioner in the form 𝑀⁻¹=𝐺^𝑇𝐺, (4.14) where 𝐺 is a sparse approximation of the inverse of the exact Cholesky factor 𝐿 of an SPD matrix 𝐴 obtained with a Frobenius norm minimization procedure. Select a lower triangular nonzero pattern 𝒮_𝐿 and find the matrix 𝐺∈𝒲_𝑆𝐿, where 𝒲_𝑆𝐿 is the set of matrices with pattern 𝒮_𝐿 such that ‖𝐼−𝐺𝐿‖_𝐹→min. (4.15) The minimization (4.15) can be performed recalling that ‖𝐼−𝐺𝐿‖_𝐹=Tr(𝐼−𝐺𝐿)(𝐼−𝐺𝐿)^𝑇=𝑛−2Tr(𝐺𝐿)+Tr𝐺𝐴𝐺^𝑇. (4.16) Deriving the right-hand side of (4.16) with respect to the nonzero entries of 𝐺 and setting to 0 yield [𝐺𝐴]_𝑖𝑗=𝑙_𝑗𝑖∀𝑖,𝑗∈𝒮_𝐿, (4.17) where the symbol [⋅]_𝑖𝑗 denotes the entry in position (𝑖,𝑗) of the matrix in the brackets. As 𝒮_𝐿 is a lower triangular pattern, the entry 𝑙_𝑗𝑖 of 𝐿 is nonzero only if 𝑖=𝑗 and the off-diagonal terms of 𝐿 are not required. Scaling (4.17) so that the right-hand side is either 1 or 0 yields 𝐺𝐴_𝑖𝑗=𝛿_𝑖𝑗∀𝑖,𝑗∈𝒮_𝐿, (4.18) where 𝐺 is the scaled factor 𝐺 and 𝛿_𝑖𝑗 the Kronecker delta. Solving (4.18) by rows leads to a sequence of dense SPD linear systems with size equal to the number of nonzeroes set for each row of 𝐺. More precisely, denoting by 𝒥 the set of column indices where the 𝑖th row 𝐠_𝑖 of 𝐺 has a nonzero entry, the following SPD systems have to be solved: 𝐴[𝒥,𝒥]𝐠_𝑖[𝒥]=𝐢_𝑚, (4.19) where i_𝑚∈ℝ^𝑚 has zero components except the 𝑚th one that is unitary, and 𝑚 is the cardinality of 𝒥. The FSAI factor 𝐺 is finally restored from 𝐺 as 𝐺=diag𝐺^−1/2𝐺 (4.20) such that the preconditioned matrix 𝐺𝐴𝐺^𝑇 has all diagonal entries equal to 1. It has been demonstrated that FSAI is optimal in the sense that it minimizes the Kaporin conditioning number of 𝐺𝐴𝐺^𝑇 for any 𝐺∈𝒲_𝑆𝐿 [131, 132]. Another remarkable property of FSAI is that the minimization in (4.15) does not require the explicit knowledge of 𝐿, which would make the algorithm unfeasible. A sketch of the procedure to build FSAI is reported in Algorithm 5. The FSAI preconditioner was later improved in several ways [133], with the aim of increasing its efficiency by dropping the smallest entries while preserving its optimal properties. Between the late 90s and early 00s further developments of FSAI were abandoned, as AINV proved generally superior. In recent years, however, a renewed interest for FSAI has started, mainly because of its high intrinsic parallel degree, for example, [134–136].

Input: Matrix 𝐴 and the nonzero pattern 𝒮𝐿 Output: Matrix 𝐺 for  𝑖=1,…,𝑛  do  Compute 𝒥={𝑗∶(𝑖,𝑗)∈𝒮𝐿}  Set 𝑚=|𝒥|  Gather 𝐴[𝒥,𝒥] from 𝐴  Solve 𝐴[𝒥,𝒥]𝐠𝑖[𝒥]=𝐢𝑚  𝐠𝑖←𝐠𝑖/√𝐠𝑖,𝑚 end

Algorithm 5 

Computation of the FSAI preconditioner.

The native FSAI algorithm has been developed for SPD matrices. A generalization to nonsymmetric matrices is possible and quite straightforward [137, 138]; however the resulting preconditioner is much less robust because the solvability of all local linear systems and the nonsingularity of the approximate inverse are no longer theoretically guaranteed. This is why with general sparse matrices SPAI and AINV are generally preferred.

4.1. Nonzero Pattern Selection and Stability Issues

A key factor for the efficiency of approximate inverse preconditioning based on the Frobenius norm minimization is the choice of the nonzero pattern 𝒮 or 𝒮_𝐿. If the nonzero pattern is not good for the problem at hand, the result is generally a failure in the Krylov solver convergence. This is due to the fact that approximate inverse techniques assume that a sparse matrix 𝑀⁻¹ can be a good approximation, in some sense, of 𝐴⁻¹. This is not obvious at all, as the inverse of a sparse matrix is structurally dense. If 𝐴 is diagonally dominant, the entries of 𝐴⁻¹ decay exponentially away from the main diagonal [139], hence a banded pattern for 𝑀⁻¹ typically provides good results. Other favorable situations include block tridiagonal or pentadiagonal matrices [140], but for a general sparse matrix choosing a good nonzero pattern is not trivial at all.

With respect to the nonzero pattern selection, approximate inverses can be defined as static or dynamic, according to whether 𝒮 and 𝒮_𝐿 are initially chosen and kept unvaried during the computation, or 𝒮 and 𝒮_𝐿 are generated by an adaptive algorithm which selects the position of the nonzero entries in an a priori unpredictable way. Typically, dynamic approximate inverses are more powerful than static ones, but their implementation especially on parallel computers can be much more complicated. The most popular and common strategy to define a static nonzero pattern for 𝑀⁻¹ is to take the nonzero pattern of a power 𝜅 of 𝐴 as is justified in terms of the Neumann series expansion of 𝐴⁻¹. In real problems, however, it is uncommon to set a 𝜅 value larger than 2 or 3, as 𝒮 can soon become too dense, so it is often preferable working with sparsified matrices. The sparsification of 𝐴^𝜅 can be done by either a prefiltration, or a postfiltration strategy, or both. For example, Chow [141] suggests dropping the entries of 𝐴 smaller than a threshold, thus obtaining the sparsified matrix 𝐴 and using the pattern of 𝐴^𝜅. An a posteriori sparsification strategy is advanced for the FSAI algorithm by Kolotilina and Yeremin [133] with the aim of reducing the cost for the preconditioner application by dropping its smallest entries. Such a dropping, however, is nontrivial and is performed so as to preserve the property that the FSAI preconditioned matrix has unitary diagonal entries, with a possibly nonnegligible computational cost. Anyway, pre- and/or postfiltration of a static approximate inverse is generally to be recommended, as the cost of the preconditioner application can be dramatically reduced with a relatively small increase of the iteration count, for example, [135, 136, 142]. By the way, this suggests that the pattern of 𝐴^𝜅 is often much larger than what would be actually required to get a good preconditioner, the difficulty relying on the a priori recognition of the most significant entries.

Dynamic approximate inverses are based on adaptive strategies which start from a simple initial guess, for example, a diagonal pattern, and progressively enlarge it until a certain criterion is satisfied. For instance, a typical dynamic algorithm is the one proposed by Grote and Huckle [117] for the SPAI preconditioner; see Algorithm 3. More specifically, they suggest to collect the indices of the nonzero components of the residual vector in a new set ℒ and to enlarge the current set 𝒥 by adding the new column indices that appear in the rows with index in ℒ. An improvement can be obtained by retaining only a few new column indices, specifically the indices 𝑗 such that the new residual norm: 𝜌_𝑗=𝐫_𝑗²₂−𝐫^𝑇_𝑗𝐴𝐞_𝑗²𝐴𝐞_𝑗²₂, (4.21) is larger than a user-specified tolerance. Then, the set ℐ is enlarged correspondingly, and a new least square problem is solved. Other dynamic techniques for generating 𝒮 have been advanced by Huckle for both nonsymmetric [143] and SPD matrices [134]. More recently, Jia and Zhu [121] have advanced a novel dynamic approximate inverse denoted as PSAI (Power Sparse Approximate Inverse), and Janna and Ferronato [144] have developed a new adaptive strategy for Block FSAI that has proved quite efficient also for the native FSAI algorithm and will be presented in more detail in the sequel.

In contrast with approximate inverses based on the Frobenius norm minimization, the AINV algorithm does not need the selection of the nonzero pattern, so it can be classified as a purely dynamic preconditioner. However, its computation turns out to be pretty similar to that of an incomplete factorization, see Algorithms 2 and 4, thus raising concerns on its numerical stability. In particular, due to incompleteness a zero pivot may result, causing a breakdown, or small pivots can trigger unstable behaviors. To avoid such difficulties a stabilized AINV preconditioner has been developed independently in [124, 145] for SPD problems. For generally nonsymmetric problems, especially with large off-diagonal entries, the numerical stability in the AINV computation may still be an issue.

Theoretically, for an arbitrary nonzero pattern 𝒮 it is not guaranteed that 𝑀⁻¹ is nonsingular. This could also happen for a few combinations of the user-specified tolerances set for generating dynamically the approximate inverse pattern. However, such an occurrence is really unlikely in practice, and typically approximate inverses prove to be very robust in problems arising from a variety of applications. According to the specific problem at hand, the choice of the most efficient type of approximate inverse may vary. For a thorough comparative study of the available options as of 1999 see the work by Benzi and Tůma [146]. Generally, SPAI is a robust algorithm able to converge also in difficult problems where incomplete factorizations fail, for example, [147], allowing the user to have a complete control of the quality of the preconditioner through the parameter 𝜀, see Algorithm 3. However, its computation is quite expensive and can be rarely competitive with AINV if only one system has to be solved so that the construction time cannot be amortized by recycling a few times the same preconditioner. For SPD problems, SPAI is not an option and FSAI is always well defined and a very stable preconditioner, even if its convergence rate may be lower than AINV. For unsymmetric matrices, however, FSAI is not reliable and may often fail. On the other hand, AINV seldom fails and generally outperforms both SPAI and FSAI.

Numerical experiments in different problems typically show that, whenever a stable ILU-type preconditioner can be computed, that is probably the most efficient preconditioner on a scalar machine [129, 148, 149]. However, approximate inverses play an important role in offering a stable and robust alternative to ILU in ill-conditioned problems. Moreover, they can become the preconditioner of choice for parallel computations, as it will be discussed in more detail in the sequence.

5. Algebraic Multigrid Methods

It is well known that the convergence rate of any Krylov subspace method preconditioned by either an incomplete factorization or an approximate inverse tends to slow down as the problem size increases. The convergence deterioration along with the increased number of operations per iteration may lead to an unacceptably high computational cost, thus limiting de facto the size of the simulated model even though large memory resources are available. This is why in recent years much work has been devoted to the so-called scalable algorithms, that is, solvers where the iteration count is insensitive to the characteristic size of the discretization mesh.

Multigrid methods can provide an answer to such a demand. Pioneering works on multigrid are due to Brandt [150, 151] who promoted these techniques for the solution of regularly discretized elliptic PDEs. The matrices arising from this kind of problem have a regular structure for which all eigenvalues and eigenvectors are known. For example, let us consider the 1D model problem: −𝑢^′′(𝑥)=𝑓(𝑥),𝑥∈[0,1], (5.1) with homogeneous Dirichlet boundary conditions discretized by centered differences over 𝑛+2 equally spaced points 𝑥_𝑖, 𝑖=0,…,𝑛+1. The grid characteristic size is therefore ℎ=1/(𝑛+1). This discretization results in the tridiagonal 𝑛×𝑛 linear system 𝐴_ℎ𝐱=𝐛, where, as is well known, the eigenvalues of 𝐴_ℎ are 𝜆_𝑘=4sin²𝑘ℎ𝜋2𝑘=1,…,𝑛, (5.2) with the associated eigenvectors 𝐰_𝑘=[sin(𝑘ℎ𝜋),sin(2𝑘ℎ𝜋),…,sin(𝑛𝑘ℎ𝜋)]^𝑇. (5.3) In other words, the 𝑖th component of 𝐰_𝑘 is actually 𝐰_𝑘,𝑖=sin𝑘𝜋𝑥_𝑖. (5.4) Let us solve the discretized model problem by a stationary method, for example, the Jacobi iteration. As the main diagonal of 𝐴_ℎ has all components equal to 2, the Jacobi iteration matrix 𝑆_ℎ is just 𝑆_ℎ𝐴=𝐼−^ℎ2 (5.5) and the 𝑘th component of the error along the direction of the 𝑘th eigenvector of 𝑆_ℎ is damped at each iteration by the reduction coefficient 𝜂_𝑘: 𝜂_𝑘=1−2sin²𝑘ℎ𝜋2. (5.6) The reduction coefficients with 𝑘 around 𝑛/2 are very small, and the related error components are damped very rapidly in a few iterations. This part of the error is often referred to as the oscillatory part corresponding to high frequency modes. However, if ℎ is small, there are also several reduction coefficients close to 1, related to both the smallest and largest values of 𝑘, that slow down significantly the stationary iteration. This part of the error is usually referred to as the smooth part corresponding to low frequency modes. Let us now consider a coarser grid where the characteristic size is doubled, that is, we simply keep the nodes of the initial discretization with an even index. Recalling (5.4) it follows that a smooth mode on the fine grid, that is, 𝑘<𝑛/4 and 𝑘>3𝑛/4, is now automatically injected into an oscillatory mode on the coarse grid. Building 𝐴_2ℎ and 𝑆_2ℎ, it is now possible to damp in a few iterations the error on the coarse grid and then project back the result on the fine grid. Similar conclusions can be drawn also using other stationary iterations, such as Gauss-Seidel, though in a much less straightforward way.

Recalling the previous observations, the basic idea of multigrid proceeds as follows. The operator 𝑆_ℎ, usually denoted as the smoother in the multigrid terminology, is applied to the fine grid discretization for 𝜈₁ iterations. The resulting fine grid residual is computed and injected into the coarse grid through a restriction operator 𝑅^𝐻_ℎ. Recalling the relationship between error and residual, the coarse grid problem is solved using the residual as known term, thus providing the coarse grid error. Such an error is then projected back to the original fine grid using a prolongation operator 𝑃^ℎ_𝐻, and the fine grid solution, obtained after the initial 𝜈₁ smoothing iterations, is corrected. Finally, the fine grid smoother is applied again until convergence. The algorithm described above is known as two-grid cycle and is the basis of multigrid techniques. Quite naturally, the two-grid cycle can be extended in a recursive way. At the coarse level, another two-grid cycle can be applied moving at a coarser grid, defining appropriate smoother, restriction, prolongation operators, and so on, until the coarsest problem is so small that it can be efficiently solved by a direct method. The recursive application of the two-grid cycle gives rise to the so-called V-cycle. Other more complex recursions can provide the so-called W-cycle; see, for instance, [152] for details.

Multigrid methods have been introduced as a solver for discretized PDEs of elliptic type, and indeed in such problems they have soon proved to be largely superior to existing algorithms. The first idea to extend their use to other applications was to look at the multigrid as a purely algebraic solver, where one has to define the smoother, restriction, and prolongation operators knowing the system matrix 𝐴 only independently of the grid and the discretized PDE it actually comes from. This strategy, known as Algebraic Multigrid (AMG), was first introduced in the mid 80s [153–155] and became soon quite popular. The second idea was to regard AMG not as a competitor with Krylov subspace solvers, rather as a potentially effective preconditioner. In fact, to work as a preconditioner it is simply sufficient to fix the number 𝜈₂ of smoothing iterations needed to reconstruct the system solution at the finest level, thus obtaining an inexact solution. The scheme of a V-cycle AMG application as a preconditioner is provided in Algorithm 6.

Input: Matrix 𝐴(𝑘), the known vector 𝐯(𝑘), the restriction and prolongation    operators 𝑅ℎ𝐻(𝑘) and 𝑃𝐻ℎ(𝑘), the smoothers 𝑆ℎ(𝑘), 𝑘, 𝜈1, 𝜈2 and 𝜇 Output: Vector 𝐮(𝑘)=𝑀−1𝐯(𝑘) if  𝑘=𝜇  then  Solve 𝐴(𝑘)𝐮(𝑘)=𝐯(𝑘) else  Apply 𝜈1 iterations of the smoother 𝑆ℎ(𝑘) to 𝐴(𝑘)𝐮(𝑘)=𝐯(𝑘)  𝐫(𝑘)=𝐯(𝑘)−𝐴(𝑘)𝐮(𝑘)  𝐫(𝑘+1)=𝑅ℎ𝐻(𝑘)𝐫(𝑘)  Call MGV(𝐴(𝑘+1),𝑅ℎ𝐻(𝑘+1),𝑃𝐻ℎ(𝑘+1),𝑆ℎ(𝑘+1),𝐫(𝑘+1),𝐞(𝑘+1))  𝐞(𝑘)=𝑃𝐻ℎ(𝑘)𝐞(𝑘+1)  𝐮(𝑘)←𝐮(𝑘)+𝐞(𝑘)  Apply 𝜈2 iterations of the smoother 𝑆ℎ(𝑘) to 𝐴(𝑘)𝐮(𝑘)=𝐯(𝑘) end

Algorithm 6 

Recursive application for 𝜇 times of the V-cycle AMG function MGV(𝐴,𝑅^𝐻_ℎ,𝑃^ℎ_𝐻,𝑆_ℎ,v,u).

Basically, AMG preconditioners vary according to the choice of both the restriction operator and the smoother, while the prolongation operator is often defined as the transposed of the restrictor. The basic idea for defining a restriction is to coarsen the native pattern of 𝐴 and prescribe an appropriate interpolation over the coarse nodes. Classical coarsening strategies have been introduced in [155, 156] which separate the variables into coarse (C-) and fine (F-) points according to the strength of each matrix connection. A point 𝑗 strongly influences the point 𝑖 if −𝑎_𝑖𝑗≥𝜃max_{𝑘=1,𝑛,𝑘≠𝑖}−𝑎_𝑖𝑘, (5.7) where 𝜃 is the user-defined strength threshold. For each point 𝑗 that strongly influences an F-point 𝑖, 𝑗 is either a C-point or is strongly influenced by another C-point 𝑘. Moreover, two C-points should not be connected to each other. An alternative successful coarsening strategy is based on the so-called smoothed aggregation (SA) technique [157], where a coarse node is defined by an aggregate of a root point 𝑖 and all the neighboring nodes 𝑗 such that 𝑎_𝑖𝑗>𝜃𝑎_𝑖𝑖𝑎_𝑗𝑗. (5.8) Another possibility is to use a variant of the independent set ordering (ISO) technique which subdivides a matrix into a 2×2 block structure and projects the native system into the Schur complement space. Using ISO or SA can lead to an elegant relationship between AMG and multilevel ILU as exploited in [98, 101, 158, 159]. More recently, a number of both parallel and scalar coarsening strategies have been also advanced, for example, [160–164]. As far as the smoother is concerned, for many years the method of choice has been the standard Gauss-Seidel iteration. With the development of parallel computers other smoothers have been introduced with the aim of increasing the parallel degree of the resulting AMG algorithm. Multicoloring approaches [165] can be used in connection to standard Gauss-Seidel, or also polynomial and C-F Jacobi smoothers [166, 167]. It is no surprise that also approximate inverses can be efficiently used as a smoother, for example, [168, 169].

The last decade has witnessed an explosion of research on AMG techniques. The key factor leading to such a great interest is basically their potential scalability with the size of the problem to be solved, in the sense that the iteration count to converge for a given problem does not depend on the number of the mesh nodes. Several theoretical results have been obtained with the aim of generalizing as much as possible AMG to nonelliptic problems, for example, [158, 170–174]. AMG techniques have been used in several different applications with promising results, such as fluid-structure interactions, meshless methods, Maxwell and Helmoltz equations, structural mechanics, sedimentary basin reconstruction, and advection-convection and reactive systems, for example, [175–185]. The above list of references is just representative and of course largely incomplete. Nonetheless, much work is still to be done in order to make AMG the method of choice. In particular, difficulties can be experienced where the system matrix arises from highly distorted and irregular meshes, or in presence of strong heterogeneities. In these cases, even the most advanced AMG strategies can fail.

6. Parallel Preconditioners

Parallel computing is widely accepted as the only pathway toward the possibility of managing millions of unknowns [186]. At the same time, hardware companies are improving the available computational resources with an increasing degree of parallelism rather than accelerating each single computing core. This is why the interest in parallel computations is continuously growing in recent years. Krylov subspace methods are in principle almost ideally parallel, as their kernels are matrix-vector and scalar products along with vector updates. Unfortunately, the same is not true for most algebraic preconditioners. There are basically two approaches for developing parallel preconditioners: the first tries to extract as much parallelism as possible from existing algorithms with the aim of transferring them to high-performance platforms, while the second is based on developing novel techniques which would not make sense on a scalar computer. Quite obviously, the former approach is easier to understand from a conceptual point of view, as the native algorithm is not modified and the difficulties are mainly technological. The latter implies an additional effort to develop new explicitly parallel methods.

ILU-based preconditioners are highly sequential in both their computation and application. Anyway, some degree of parallelism can be achieved through graph coloring techniques. These strategies have been known from a long time by numerical analysts, for example, [6], and used in the context of relaxation methods. The aim is to build a partition of the matrix graph such that all vertices in the same group form an independent set, so that all unknowns in the same subset can be solved in parallel. Other useful orderings are based on domain decomposition strategies minimizing the edges connecting different subdomains, with a local ordering applied also to each unknown subset. Parallel ILU implementations based on these concepts have been developed in [187–190]; however it was soon made clear that generally such algorithms could not have a promising scalability.

In contrast, approximate inverses are intrinsically much more parallel than ILU factorizations, as they can be applied by a matrix-vector product instead of forward and backward substitutions. The construction of an approximate inverse, however, can be difficult to parallelize. For example, the AINV Algorithm 4 proceeds by columns and is conceptually pretty similar to an ILU factorization. A parallel AINV construction can be obtained by means of graph partitioning techniques which decompose the adjacency graph of 𝐴 in a number of subsets of roughly the same size with a minimal number of edge cuts [191]. In such a way the inner variables of each subset can be managed independently by each processor, while communications occur only with the boundary variables. By distinction, the construction of an approximate inverse based on a Frobenius norm minimization is naturally parallel if a static nonzero pattern is defined a priori. For example, the computation of both static SPAI and FSAI consists of a number of independent least-square problems and dense linear systems, respectively, which can be trivially partitioned among the available processors. Hence, static SPAI and FSAI are almost perfectly parallel algorithms. However, when the nonzero pattern is defined dynamically parallelization is no longer straightforward, as the operation of gathering the required submatrix of 𝐴, see Algorithms 3 and 5, gives rise to unpredictable communications. For example, details of the nontrivial parallel implementation of the dynamic SPAI algorithm are provided in [147].

The theoretical scalability properties of AMG make it very attractive for parallel computations. This is why in the last few years the research on AMG techniques has concentrated on high-performance massively parallel implementations. The main difficulties may arise in parallelizing the Gauss-Seidel smoother and the coarsening stage. As mentioned before, these problems can be overcome by using naturally parallel smoothers, such as Jacobi, relaxed Jacobi, polynomial or a static approximate inverse smoother, for example, [162, 166], and parallel coarsening strategies, for example, [160, 162, 192, 193]. Currently, efficient parallel implementations of AMG are already available in several software packages, for example, Hypre [194] and Trilinos [195].

An alternative strategy for developing parallel preconditioners relies on building new algorithms which should consist of matrix-vector products and local triangular solves only. Perhaps the earliest technique of this kind belongs to the class of the polynomial preconditioners which was first introduced by Cesari as back as in 1937 [4], even though obviously not in the context of Krylov subspace methods and supercomputing. Modern polynomial preconditioners have been developed in [196–198]. A polynomial preconditioner is defined as 𝑀⁻¹=𝑠(𝐴), (6.1) where 𝑠∈ℙ_𝑘 is a polynomial of degree not exceeding 𝑘. The preconditioned system becomes: 𝑠(𝐴)𝐴𝐱=𝑠(𝐴)𝐛, (6.2) where 𝑠(𝐴) and 𝐴 commute, so that right and left preconditionings coincide. Quite obviously, 𝑠(𝐴)𝐴 is never formed explicitly and its application to a vector is simply performed by a sequence of matrix-vector products. Therefore, polynomial preconditioners can be ideally implemented in a parallel environment. There are different polynomial preconditioners according to the choice of 𝑠(𝐴). One of the simplest options, advanced in [196], is using the Neumann expansion 𝐼+𝑁+𝑁²+⋯+𝑁^𝑙, where 𝑁=𝐼−𝜔𝐴 (6.3) and 𝜔 is a scaling parameter. The degree 𝑙≤𝑘 cannot be too large, and this can limit the effectiveness of this preconditioner. Another option, advanced in [197], is to select a somewhat “optimal” polynomial, where optimality is evaluated in some sense. A natural choice is prescribing that the eigenspectrum of the preconditioned matrix is as close as possible to that of the identity, that is, finding 𝑠∈ℙ_𝑘 such that max_{𝜆∈𝜎(𝐴)}|1−𝜆𝑠(𝜆)|→min, (6.4) where 𝜎(𝐴) is the eigenspectrum of 𝐴. The min-max problem above can be addressed by using Chebyshev polynomials but cannot be solved exactly, as 𝜎(𝐴) is generally not known. An approximate minimization can be done replacing 𝜎(𝐴) with a continuous set which encloses 𝜎(𝐴). A third option is to compute 𝑠 as the polynomial that minimizes ‖1−𝜆𝑠(𝜆)‖_𝑤, (6.5) where the symbol ‖⋅‖_𝑤 denotes the 2-norm induced by the inner product on the space ℙ_𝑘 with respect to some nonnegative weight function defined in the interval [𝛼,𝛽]: 𝑝,𝑞=^𝛽_𝛼𝑝(𝜆)𝑞(𝜆)𝑤(𝜆)𝑑𝜆. (6.6) The resulting polynomial is known as least-squares and was actually introduced by Stiefel in the solution of eigenproblems [199]. According to Saad [198], least-squares polynomials tend to perform better than Chebyshev and Neumann. However, a major drawback for the use of this kind of preconditioners stems from the need at least for estimates of the smallest and largest eigenvalues of 𝐴 and often their performance is quite poor in terms of convergence rate [200].

Another popular parallel preconditioner is based on the Additive Schwarz (AS) methods [201]. The idea is to use a domain decomposition approach for dividing the variables into a number of possibly overlapping subdomains and to address separately each subdomain at a local level (Figure 3). On the overlaps the contribution of each subdomain is added or treated in some special way by averaging techniques. AS preconditioners can be easily parallelized by assigning each subdomain to a single processor. For instance, a sequential ILU factorization can be computed for each subdomain, thus miming a kind of parallel ILU. If no overlaps are allowed for, AS preconditioners with inner ILU factorizations simply reduce to an incomplete Block Jacobi. This class of preconditioners has experienced some success in the 90s as it is conceptually perfectly parallel. Unfortunately, the quality of the preconditioner deteriorates as the number of processors, that is, subdomains, grows, because the size of each subset progressively decreases. This causes an increase of the iteration count which can also completely offset the advantage of using a larger number of processors.

Figure 3 

Matrix partitioning into possibly overlapping subdomains.

A novel recent parallel preconditioner which tries to combine the positive features of both approximate inverses based on the Frobenius norm minimization and AS methods is Block FSAI [202]. The basic idea exploits the fact that in the Frobenius norm minimization the preconditioned matrix can be actually forced to resemble any target matrix 𝑇, by computing a pseudoapproximate inverse 𝑀⁻¹ such that 𝑇−𝐴𝑀⁻¹_𝐹→min (6.7) for all matrices having a prescribed nonzero pattern 𝒮. This concept, originally introduced in [120], can be used to choose a target matrix 𝑇 that is particularly attractive for parallel computing. The idea of Janna et al. [202] is to generalize the FSAI approach using as a target a block diagonal matrix. As FSAI is a robust preconditioner for SPD problems, Block FSAI has been originally developed for this class of matrices, although a generalization to nonsymmetric indefinite linear systems has also been attempted [203]. So let 𝐴 be SPD with 𝒮_𝐿 and 𝒮_𝐵𝐷 a sparse lower triangular and a dense block diagonal nonzero 𝑛×𝑛 pattern, respectively. Denote by 𝑛_𝑏 the number of diagonal blocks and 𝑚_𝑖𝑏 the size of the 𝑖_𝑏th block, and let 𝐷 be an arbitrary full-rank matrix with nonzero pattern 𝒮_𝐵𝐷. The Block FSAI preconditioner of 𝐴 is defined as the product 𝐹^𝑇𝐹, where 𝐹 is a lower block triangular factor with nonzero pattern 𝒮_𝐵𝐿=𝒮_𝐵𝐷∪𝒮_𝐿 such that ‖𝐷−𝐹𝐿‖_𝐹→min (6.8) with 𝐿 the exact lower Cholesky factor of 𝐴. As 𝐷 is arbitrary, it goes without saying that 𝐹 as defined in (6.8) is also. Similarly to the FSAI algorithm, minimization of (6.8) yields a linear relationship for the 𝑖th row 𝐟_𝑖 of 𝐹: 𝐴[𝒥,𝒥]𝐟_𝑖[𝒥]=𝐯, (6.9) where 𝒥 is the set of column indices in the 𝑖th row of 𝐹 with nonzero entries and 𝐯 is the null vector except for the last 𝑚_𝑖𝑏 components which are arbitrary, being 𝑖_𝑏 the block index the row 𝑖 belongs to (Figure 4). For the system (6.9) to have a unique solution 𝑚_𝑖𝑏 components of 𝐟_𝑖[𝒥] can be set arbitrarily because of the arbitrariness of 𝐷. According to [202] the most convenient choice is to prescribe all the components of 𝐟_𝑖[𝒥] falling within the 𝑖_𝑏th diagonal block to be null, with the exception of the diagonal entries set to 1. This implies that 𝐹 is actually a unit lower triangular matrix with structure 𝐹𝐹=𝐼0⋯0₂₁𝐹𝐼⋮⋮⋱⋮_𝑛𝑏1𝐹_𝑛𝑏2⋯𝐼. (6.10) Similarly to the static FSAI algorithm, the factor 𝐹 can be efficiently built in parallel as all systems (6.9) are independent. According to (6.8), the preconditioned matrix 𝐹𝐴𝐹^𝑇 should resemble 𝐷𝐷^𝑇, that is, a block diagonal matrix for any arbitrary 𝐷. In other words, 𝐹𝐴𝐹^𝑇 has the block structure 𝐹𝐴𝐹^𝑇=𝐵₁𝑅_1,2⋯𝑅_1𝑛𝑏𝑅^𝑇₁₂𝐵₂⋯𝑅_2𝑛𝑏𝑅⋮⋱⋮_{𝑇1𝑛𝑏}𝑅^𝑇_2𝑛𝑏⋯𝐵_𝑛𝑏, (6.11) where 𝑅_𝑖𝑗 are residual blocks and 𝐵_𝑖𝑏 tend to the diagonal blocks of 𝐷𝐷^𝑇. As 𝐷 is arbitrary, there is no reason for 𝐹𝐴𝐹^𝑇 to be better than 𝐴 in a CG iteration, so it is necessary to precondition 𝐹𝐴𝐹^𝑇 again. Assuming that the residual off-diagonal blocks in (6.11) are close to 0, an effective choice could be using as “outer” preconditioner a nonoverlapped AS method containing a local preconditioner for each block 𝐵_𝑖𝑏. Writing the “outer” preconditioner in the factored block diagonal form 𝐽^−𝑇𝐽⁻¹, the resulting preconditioned matrix reads 𝐽⁻¹𝐹𝐴𝐹^𝑇𝐽^−𝑇=𝑊𝐴𝑊^𝑇 (6.12) with the final Block FSAI preconditioner: 𝑀⁻¹=𝑊^𝑇𝑊=𝐹^𝑇𝐽^−𝑇𝐽⁻¹𝐹. (6.13)

(a)

(b)

Figure 4 

Elements belonging to 𝐟𝑖[𝒥](a) and 𝐴[𝒥,𝒥](b).

The Block FSAI preconditioner described above can be improved by defining 𝒮_𝐵𝐿 dynamically. An adaptive algorithm for the 𝒮_𝐵𝐿 pattern identification can be implemented starting from the theoretical optimal properties of Block FSAI. Janna and Ferronato [144] have demonstrated that under the hypothesis that the 𝐽^−𝑇𝐽⁻¹ contains the exact inverse of each diagonal block 𝐵_𝑖𝑏 the Kaporin conditioning number 𝛽 of the preconditioned matrix 𝑊𝐴𝑊^𝑇 [132]: 𝛽𝑊𝐴𝑊^𝑇=tr𝑊𝐴𝑊^𝑇𝑛det𝑊𝐴𝑊^𝑇1/𝑛, (6.14) satisfies the following bound: 1≤𝛽𝑊𝐴𝑊^𝑇≤𝐶𝜓(𝐹), (6.15) where 𝐶 is a constant scalar independent of 𝑊 and 𝜓(𝐹)=^𝑛_𝑖=1𝐹𝐴𝐹^𝑇_𝑖𝑖^1/𝑛. (6.16)

Block FSAI has the theoretical property of minimizing 𝜓(𝐹) for any given pattern 𝒮_𝐵𝐿. The basic idea of the adaptive algorithm is to select the off-block diagonal nonzero entries in any row 𝐟_𝑖 of 𝐹 so as to reduce 𝜓(𝐹) as much as possible. This is feasible because each factor [𝐹𝐴𝐹^𝑇]_𝑖𝑖 turns out to be a quadratic form depending on 𝐟_𝑖 only. Denoting by 𝐟_𝑖̃ the subvector of 𝐟_𝑖 including the off-block diagonal entries only, 𝐴_𝑖𝑏 the square submatrix of 𝐴 built from the first to the 𝑚th row/column, with 𝑚 the sum of size of the first (𝑖_𝑏−1) diagonal blocks of 𝒮_𝐵𝐷, and 𝐚_𝑖 the subrow of 𝐴 with the first 𝑚 elements of the 𝑖th row (Figure 5), each factor [𝐹𝐴𝐹^𝑇]_𝑖𝑖 in (6.16) reads 𝐹𝐴𝐹^𝑇_𝑖𝑖=𝐟^𝑇_𝑖̃𝐴_𝑖𝑏𝐟_𝑖̃+2𝐟^𝑇_𝑖̃𝐚_𝑖+𝑎_𝑖𝑖. (6.17) Minimizing every [𝐹𝐴𝐹^𝑇]_𝑖𝑖, 𝑖=1,…,𝑛, is equivalent to minimizing 𝜓(𝐹). The adaptive pattern search for 𝐹 can be therefore implemented as follows. Start from an initial guess 𝒮⁰_𝐹 for the nonzero pattern of the off-block diagonal part of 𝐹, for example, the empty pattern, and compute the gradient 𝐠_𝑖 of each quadratic form [𝐹𝐴𝐹^𝑇]_𝑖𝑖: 𝐠_𝑖=2𝐴_𝑖𝑏𝐟_𝑖̃+𝐚_𝑖. (6.18) Then, add to 𝒮⁰_𝐹 the position 𝑗 into each row 𝑖 corresponding to the largest component of 𝐠_𝑖 computed in (6.18), so as to obtain an augmented pattern 𝒮¹_𝐹. After computing the new factor 𝐹 over 𝒮¹_𝐹, the procedure above can be iterated in order to build 𝒮²_𝐹, 𝒮³_𝐹, and so on. The search into each row can be stopped when either a maximum number 𝑘_max of entries are added to the initial pattern or the relative variation of [𝐹𝐴𝐹^𝑇]_𝑖𝑖 in two consecutive steps is smaller than a prescribed tolerance 𝜖_𝐹. The construction of the dynamic Block FSAI is summarized in Algorithm 7. The Block FSAI computed with the dynamic algorithm described above is known as Adaptive Block FSAI (ABF) and has proved very effective even if the outer preconditioner does not contain the exact inverse of the diagonal blocks 𝐵_𝑖𝑏. As any dynamic algorithm, the parallel construction of ABF is not straightforward. However, an almost ideal efficiency can be obtained in massively parallel computations as well. For more details, see [144].

Input: Matrix 𝐴, the initial pattern 𝒮0𝐹 Output: Matrix 𝐹 for  𝑖=1,…,𝑛  do  Extract 𝒥0 from 𝒮0𝐹  if  𝒥0≠∅  then   Gather 𝐴[𝒥0,𝒥0] and 𝐴[𝒥0,𝑖] from 𝐴   Solve 𝐴[𝒥0,𝒥0]𝐟𝑖[𝒥0]=−𝐴[𝒥0,𝑖]  end  Set 𝑘=0 and 𝜖𝑖>𝜖𝐹  while  (𝑘≤𝑘max) and (𝜖𝑖<𝜖𝐹)  do   𝑘←𝑘+1   Compute 𝐠𝑖=2(𝐴𝑖𝑏𝐟𝑖̃+𝐚𝑖)   Update 𝒥𝑘−1 by adding the index of the largest component of 𝐠𝑖   Gather 𝐴[𝒥𝑘,𝒥𝑘] and 𝐴[𝒥𝑘,𝑖] from 𝐴   Solve 𝐴[𝒥𝑘,𝒥𝑘]𝐟𝑖[𝒥𝑘]=−𝐴[𝒥𝑘,𝑖]   Compute the diagonal term [𝐹𝐴𝐹𝑇]𝑘𝑖𝑖 at the current step   Compute 𝜖𝑖=([𝐹𝐴𝐹𝑇]𝑘𝑖𝑖−[𝐹𝐴𝐹𝑇]𝑘−1𝑖𝑖)/[𝐹𝐴𝐹𝑇]𝑘−1𝑖𝑖  end end

Algorithm 7 

Dynamic construction of the Block FSAI preconditioner.

(a)

(b)

Figure 5 

Matrix 𝐴𝑖𝑏 and vectors 𝐟𝑖̃ and 𝐚𝑖 used in the [𝐹𝐴𝐹𝑇]𝑖𝑖 and 𝐠𝑖 computation.

Block FSAI coupled with a block diagonal IC factorization, for example, as implemented in [74, 84, 86, 90], has proved to be a very efficient parallel preconditioner for both SPD linear systems [144, 202] and eigenproblems [204, 205] arising from different areas, such as fluid dynamics, structural mechanics, and contact problems. The combination of FSAI with block diagonal IC joins the parallelization degree of the former with the effectiveness of the latter. Similarly to standard incomplete Block Jacobi algorithms, the iteration count tends to grow increasing the number of computing cores. However, the convergence deterioration is much smoothed by the Block FSAI effect which always outperforms the standard FSAI algorithm. To improve the preconditioner scalability for massively parallel simulations the use of graph partitioning and/or domain decomposition strategies can be of great help [206, 207].

7. Saddle-Point Problems

In recent years the interest in the solution of saddle-point problems has grown, with many contributions from different fields. The main reason is the great variety of applications requiring the solution of linear systems with a saddle-point matrix. For example, such problems arise in the discretization of compressible and incompressible Navier-Stokes equations in computational fluid dynamics [208, 209], constrained optimization [210, 211], electromagnetism [212, 213], mixed FE approximations in fluid and solid mechanics [208, 214], and coupled consolidation problems in geomechanics [58].

A saddle-point problem is defined as a linear system where the matrix 𝐴 has the structure 𝐴=𝐾𝐵₁𝐵₂−𝐶, (7.1) where 𝐾∈ℝ^𝑛1×𝑛1, 𝐵₁∈ℝ^𝑛1×𝑛2, 𝐵₂∈ℝ^𝑛2×𝑛1, 𝐶∈ℝ^𝑛2×𝑛2, and 𝐾, 𝐵₁, and 𝐵₂ are nonzero. Typically, but not necessarily, 𝑛₁≥𝑛₂. Moreover, 𝐶 is symmetric and positive semidefinite, and often 𝐶=0, while in most applications 𝐾 is SPD and 𝐵₂=𝐵^𝑇₁, that is, 𝐴 is symmetric indefinite. For example, a saddle-point matrix 𝐴 arises naturally in an equality-constrained quadratic problem: 𝐽(𝐱)=12𝐱^𝑇𝐴𝐱−𝐟^𝑇𝐱→minsubjectto𝐵𝐱=𝐠 (7.2) with 𝐶=0 and 𝐵₂=𝐵^𝑇₁=𝐵, solved with the aid of Lagrangian multipliers. A nonzero 𝐶 may arise in some interior point methods [210, 211], compressible Navier-Stokes equations [215], or coupled consolidation problems [65]. Frequently the norm of 𝐶 is much smaller than that of 𝐾. Because of the great variety of applications leading to saddle-point matrices, it is no surprise that several different solution methods have been developed. For an exhaustive, review, see Benzi et al. [216]. In the context of the present paper, the attention will be obviously restricted to preconditioned iterative methods.

A natural way to solve the problem 𝐾𝐵^𝑇𝐱𝐵−𝐶𝐱¹₂=𝐛₁𝐛₂ (7.3) is to reduce the system (7.3) to the Schur complement form. Solve for 𝐱₁ in the upper set of equations: 𝐱₁=𝐾⁻¹𝐛₁−𝐵^𝑇𝐱₂, (7.4) and substitute (7.4) in the lower set: 𝐵𝐾⁻¹𝐵^𝑇+𝐶𝐱₂=𝐵𝐾⁻¹𝐛₁−𝐛₂. (7.5) The matrix 𝑆=𝐵𝐾⁻¹𝐵^𝑇+𝐶 at the left-hand side of (7.5) is the Schur complement of system (7.3). Hence, the global solution 𝐱 can be obtained by solving the Schur complement system (7.5) and substituting the result back in (7.4) which implies solving another system with the matrix 𝐾. This pair of SPD linear systems can be elegantly solved by partitioned schemes, as those proposed in [217, 218] in the context of coupled diffusion and consolidation problems, resulting in a double CG iteration where two distinct preconditioners for 𝐾 and 𝑆 are needed. A usually more efficient alternative is to solve the systems (7.4) and (7.5) inexactly and use the procedure described above as a preconditioner in a Krylov subspace method. This is what is traditionally referred to as constraint preconditioning because the idea was first introduced in the solution of equality-constrained optimization problems such as (7.2) [219, 220].

Quite obviously, there are several different constraint preconditioners according to how 𝐾⁻¹ and 𝑆⁻¹ in (7.4) and (7.5) are approximated. The most effective choice turns out to be strongly dependent on the specific problem at hand. A popular strategy consists of replacing 𝐾⁻¹ with 𝑀_𝐾⁻¹ and solving exactly the Schur complement system (7.5) [219–222]. Such an approach can be convenient if 𝑛₂ is much smaller than 𝑛₁, so that the system (7.5) can be efficiently solved by a direct method. Otherwise, an inner CG iteration can be set up exiting the procedure even at a relatively large residual [130]. In this case the resulting preconditioner reads 𝑀⁻¹_ECP=𝑀_𝐾𝐵^𝑇𝐵−𝐶⁻¹ (7.6) and is usually denoted as Exact Constraint Preconditioner (ECP). ECP has several interesting theoretical properties which have been studied in detail by several authors. For instance, an exhaustive eigenanalysis of the preconditioned matrix can be found in [223]. In particular, denoting with 𝛼_𝐾 and 𝛽_𝐾 the smallest and the largest eigenvalue, respectively, of 𝑀_𝐾⁻¹𝐾, Durazzi and Ruggiero [224] have shown that the eigenvalues 𝜆 of 𝑀⁻¹_ECP𝐴 are either one, with multiplicity at least 𝑛₂, or real positive and bounded by 𝛼_𝐾≤𝜆≤𝛽_𝐾. Quite interestingly, this result can lead the classical CG algorithm to converge even with the indefinite matrix (7.3), provided that the last 𝑛₂ components of the initial residual are zero. This can be simply accomplished by choosing as initial guess 𝑀⁻¹_ECP𝐛.

Despite its remarkable properties, ECP has at least two severe limitations. First, 𝑀_𝐾⁻¹ must be known explicitly to allow for the computation of 𝑆 and must be very sparse in order to preserve a workable sparsity for 𝑆 as well. Second, although an accurate solution of (7.5) is not really needed, solving it at each preconditioner application can represent a significant computational burden for the overall scheme. This is why ECP is often used with 𝑀_𝐾⁻¹=[diag(𝐾)]⁻¹, which can be quite detrimental for the solver convergence in case 𝐾 is not diagonally dominant or is ill-conditioned, such as in coupled consolidation problems [130].

To make the preconditioner application cheaper, a substitute for (7.5) can be solved replacing 𝑆 with another approximation 𝑀_𝑆 as easy to invert as possible. Moreover, in some cases it may also be possible to approximate 𝑆 without computing it explicitly, thus further reducing the constraint preconditioner cost. Obviously the theoretical properties of ECP no longer hold true. Recalling the factorization (2.12) by Sylvester's theorem of inertia, the new preconditioner is 𝑀⁻¹_ICP=𝐼−𝑀_𝐾⁻¹𝐵^𝑇𝑀0𝐼^𝐾−100−𝑀_𝑆⁻¹𝐼0−𝐵𝑀_𝐾⁻¹𝐼=𝑀_𝐾𝐵^𝑇𝐵𝐵𝑀_𝐾⁻¹𝐵^𝑇−𝑀_𝑆⁻¹ (7.7) and is denoted as Inexact Constraint Preconditioner (ICP). This class of preconditioners has been introduced in [130, 225, 226] using [diag(𝐾)]⁻¹ or AINV for 𝑀_𝐾⁻¹ and a zero fill-in IC factorization for 𝑀_𝑆⁻¹. The choice for 𝑀_𝑆⁻¹ is justified by the fact that usually 𝑆 is well conditioned. Moreover, this allows for splitting the central block diagonal factor in (7.7) into the product of two matrices, thus providing a factorized form for 𝑀⁻¹_ICP. A block triangular variant of ICP can be easily deduced from (7.7) by neglecting either the upper or the lower outer factors: 𝑀⁻¹_TICP=𝐼−𝑀_𝐾⁻¹𝐵^𝑇𝑀0𝐼^𝐾−100−𝑀_𝑆⁻¹=𝑀_𝑘𝐵^𝑇0−𝑀_𝑆⁻¹. (7.8) The theoretical properties of block triangular ICP have been investigated by Simoncini [227] and effectively used in the solution of the discretized Navier-Stokes equations [228–230].

Similar to ECP, ICP as well requires an explicit approximation 𝑀_𝐾⁻¹ to build the Schur complement. This is still a limitation in some problems, such as coupled consolidation [130], where even AINV can be a poor approximation of 𝐾⁻¹ leading to an unacceptably high iteration count to converge. Actually, this is only partially connected with the quality of AINV itself, rather it is the need for the explicit construction of the Schur complement matrix that calls for a reduced fill-in of 𝑀_𝐾⁻¹. In order to remove this inconvenience an implicit approximation of 𝐾⁻¹ can be used, such as a stabilized IC factorization with partial fill-in: 𝑀_𝐾⁻¹=𝐿_𝐾^−𝑇𝐿_𝐾⁻¹. (7.9) The new Schur complement reads 𝑆=𝐵𝐿_𝐾^−𝑇𝐿_𝐾⁻¹𝐵^𝑇+𝐶 (7.10) but cannot be computed explicitly. Although the 𝑆 implicit form (7.10) can still be used to perform a product between 𝑆 and a vector, and consequently one may think of solving the related inner system by an iterative method, a suitable preconditioner for 𝑆 is not easily available. Alternatively, an explicit approximation of 𝑆 can be computed, for example, using the AINV of 𝐾⁻¹, namely 𝑆=𝐵𝑍𝑍^𝑇𝐵^𝑇+𝐶, (7.11) where 𝑍 is the upper triangular factor resulting from the AINV Algorithm 4, and then an IC factorization of 𝑆 is performed: 𝑀_𝑆⁻¹=𝐿_𝑆^−𝑇𝐿_𝑆⁻¹. (7.12) In this way the quality of 𝑀_𝐾⁻¹ is improved, but generally an additional approximation in 𝑀_𝑆⁻¹ is introduced. However, in several problems, especially those arising from coupled consolidation where such a constraint preconditioning variant has been introduced [148, 231], the overall algorithm has a beneficial effect. This ICP variant blending two different approximations for 𝐾⁻¹ in the same scheme reads 𝑀⁻¹_MCP=𝐼−𝐿_𝐾^−𝑇𝐿_𝐾⁻¹𝐵^𝑇𝐿0𝐼^𝐾−𝑇𝐿_𝐾⁻¹00−𝐿_𝑆^−𝑇𝐿_𝑆⁻¹𝐿𝐼0−𝐵_𝐾^−𝑇𝐿_𝐾⁻¹𝐼=𝐿_𝐾^−𝑇−𝐿_𝐾^−𝑇𝐺^𝑇0𝐿_𝑆^−𝑇𝐿_𝐾⁻¹0𝐺𝐿_𝐾⁻¹−𝐿_𝑆⁻¹, (7.13) where 𝐺=𝐿_𝑆⁻¹𝐵𝐿_𝐾^−𝑇, and is denoted as Mixed Constraint Preconditioner (MCP). Recently an improved MCP version has been advanced by properly relaxing the Schur complement approximation [61].

It is quite well accepted that constraint preconditioners outperform any other “global” preconditioning approach applied to saddle-point matrices. This remarkable property has attracted the theoretical interest of several researchers who have investigated the spectral properties of constraint-preconditioned saddle-point matrices [225, 227, 231–233], obtaining quite accurate bounds for the eigenvalues of 𝑀⁻¹𝐴. For instance, in the case of MCP in (7.13) Ferronato et al. [231] have shown that the eigenvalues 𝜆 of 𝑀⁻¹_MCP satisfy the following inequality: |𝜆−1|≤𝑒_𝑘1+𝑔²+𝑒_𝑠, (7.14) where 𝑒_𝑘=𝐼−𝐿_𝐾⁻¹𝐾𝐿_𝐾^−𝑇,𝑒_𝑠=𝐼−𝐿_𝑆⁻¹𝑆𝐿_𝑆^−𝑇,𝑔=‖𝐺‖. (7.15) In other words, 𝑒_𝑘 and 𝑒_𝑠 are a measure of the error made when using the approximations 𝑀_𝐾⁻¹ and 𝑀_𝑆⁻¹ in place of 𝐾⁻¹ and 𝑆⁻¹, respectively, while 𝑔 provides an estimate of the strength of coupling between the upper and lower set of equations in (7.3). Inequality (7.14) provides a useful and relatively cheap indication on how the quality of 𝑀_𝐾⁻¹ and 𝑀_𝑆⁻¹ affects the overall convergence rate of MCP. In particular, note that the quality of 𝑀_𝐾⁻¹ is generally more important than that of 𝑀_𝑆⁻¹, with the difference becoming larger as the coupling becomes stronger.

Quite naturally, there is the distinct interest to implement constraint preconditioners on parallel computers. The main difficulties for developing efficient parallel variants of constraint preconditioners are twofold. First, the constraint preconditioning concept is inherently sequential, as it involves the preliminary computation of the Schur complement and then the subsequent approximate solution to two inner systems. Thus, a standard parallel approach where the system matrix is distributed among the processors as stripes of contiguous rows is not feasible, because all the processors owning the second block of unknowns are idle when the other processors perform the operation on the first block, and viceversa. Second, efficient parallel approximations of 𝐾⁻¹ and 𝑆⁻¹ have to replace incomplete factorizations. For example, in incompressible Navier-Stokes problems parallel AMG techniques have been efficiently used to approximate the 𝐾 block which arises from a discretized Laplace operator [234]. In coupled consolidation better results have been obtained by using a double FSAI [235], parallel multilevel variants [236], and Block FSAI [237], which appears to be currently the most promising approach in the context of geomechanical problems.

8. Conclusions and Future Perspectives

The development and implementation of efficient preconditioners are the key factors for improving the performance and robustness of Krylov subspace methods in the solution of large sparse linear systems. Such as issue is a central task in a great number of different applications, not all necessarily related to PDEs, so it is no surprise that much research has been devoted to it. After a time where most efforts were focused on direct and iterative solvers, the last two decades can be denoted as the “preconditioning era.” The number of papers, projects, and international conferences addressing this topic has largely exceeded those aimed at the development of new solvers, and this trend is not likely to fade in the next few years. This is mainly due to the fact that an optimal general-purpose preconditioner does not exist and that any specific problem can afford the opportunity to develop an ad hoc strategy.

The present paper is intended to offer an overview of the most important algebraic preconditioners available today for the solution of sparse linear systems. Hence, it cannot exhaust the subject of preconditioning as the entire class of problem-specific algorithms has been omitted. Algebraic preconditioners have been chosen for their “universality,” in that they can be effectively used as black-box tools in general-purpose codes requiring the knowledge of the system matrix only. A classification of algebraic preconditioners has been attempted based on the most significant representatives of each group along with their prominent features: (i)incomplete factorizations have been the earliest algebraic preconditioners massively used to accelerate Krylov subspace solvers. The several weak points related to robustness and numerical stability have stimulated the development of many different variants meant to both improve their performance and extend the application area. Although some advances are still under way, this class of preconditioners appears to be quite mature with the related strong and weak points very well known. Incomplete factorizations can be very effective and in many cases are still the most performant preconditioners, but their intrinsic lack of robustness and especially the very reduced parallel degree are probably going to bound their role in the next future; (ii)approximate inverses were originally developed with the aim of resolving the main difficulties of ILU-based preconditioners. Between the late 90s and early 00s they looked as a very promising approach especially in view of their robustness and potentially high parallel degree. Indeed, these preconditioners have proved very robust and so well suited for black-box solvers and in some cases almost perfectly parallel. In a sequential run, however, their performance can be lower than that of an ILU-based preconditioner. At present, the class of approximate inverses appears to be quite a mature field of research, with the most significant breakthroughs advanced more than 10 years ago. Nonetheless, their influence is probably going to be still very important in the next few years; (iii)the potentially high scalability of algebraic multigrid techniques has promoted much research in this area. In contrast with both incomplete factorizations and approximate inverses, AMG can allow for a stable iteration count to converge independently of the characteristic mesh size, that is, the overall system dimension, thus making AMG a potentially ideally scalable tool. Currently, AMG appears as one of the more active research areas, with possible contaminations from incomplete factorizations and approximate inverses in the choice of the smoother. In several applications arising from discretized PDEs with severe mesh distortions, however, AMG can still lack robustness, and more research is still needed to make it the algorithm of choice in the next future; (iv)the current hardware development is leading to a much more pronounced parallel degree of any computer. Besides the parallelization of existing algorithms, this trend is going to enforce the development of totally new techniques which can make sense in a parallel computing environment only. This is why parallel preconditioners are now emerging as a class of algorithms which are likely to expand in the next future. These tools can take advantage of important contributions from all the previous classes and can also relaunch some techniques, such as polynomial preconditioners or FSAI-based approximate inverses, which had been somewhat forgotten over the last years.

Acknowledgments

The author is very thankful to Giuseppe Gambolati and Carlo Janna for their careful reading of the paper and helpful comments.