- Home
- Documents
*Preconditioning for Sparse Linear Systems at the Dawn of ... solvers of choice were obviously direct*

prev

next

out of 49

View

0Download

0

Embed Size (px)

International Scholarly Research Network ISRN Applied Mathematics Volume 2012, Article ID 127647, 49 pages doi:10.5402/2012/127647

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

Massimiliano Ferronato

Dipartimento ICEA, Università di Padova, Via Trieste 63, 35121 Padova, Italy

Correspondence should be addressed to Massimiliano Ferronato, ferronat@dmsa.unipd.it

Received 1 October 2012; Accepted 22 October 2012

Academic Editors: J. R. Fernandez, Q. Song, S. Sture, and Q. Zhang

Copyright q 2012 Massimiliano Ferronato. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

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 algo- rithms requiring the knowledge of the systemmatrix only, independently of the specific problem it arises from. Along with the traditional distinction between incomplete factorizations and approxi- mate 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 “precondi- tioner” is the mathematical operator responsible for such a transformation. In the context of linear systems of equations,

Ax � b, �1.1�

where A ∈ Rn×n and x,b ∈ Rn, the preconditioner is a matrix which transforms �1.1� into an equivalent problem for which the convergence of an iterative method is much faster.

2 ISRN Applied Mathematics

Preconditioners are generally used when the matrix A 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 math- ematician 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 cur- rently 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 enter- tainment, as one can think about many other things while carrying out the repetitive opera- tions 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 likeWorldWar 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 computingmachines, which later led to themodern 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 systemmatrix 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 m steps, if m is the number

ISRN Applied Mathematics 3

of distinct eigenvalues of the system matrix A. 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 m. 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 newmethod 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 inmany 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 prob- lems 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 n steps, being n the size of A. 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 Kryl