Templates for the Solution of Linear Systems: Building Blocks for Iterative Methods

The front matter includes the title page, copyright page, author's affiliations, how to use this book, acknowledgments, TOC, list of symbols, and list of figures.

Which of the following statements is true?

• Users want “black box” software that they can use with complete confidence for general problem classes without having to understand the fine algorithmic details.

• Users want to be able to tune data structures for a particular application, even if the software is not as reliable as that provided for general methods.

It turns out both are true, for different groups of users.

Traditionally, users have asked for and been provided with black box software in the form of mathematical libraries such as LAPACK, LINPACK, NAG, and IMSL. More recently, the high-performance community has discovered that they must write custom software for their problem. Their reasons include inadequate functionality of existing software libraries, data structures that are not natural or convenient for a particular problem, and overly general software that sacrifices too much performance when applied to a special case of interest.

The term “iterative method” refers to a wide range of techniques that use successive approximations to obtain more accurate solutions to a linear system at each step. In this book we will cover two types of iterative methods. Stationary methods are older, simpler to understand and implement, but usually not as effective. Nonstationary methods are a relatively recent development; their analysis is usually harder to understand, but they can be highly effective. The nonstationary methods we present are based on the idea of sequences of orthogonal vectors. (An exception is the Chebyshev iteration method, which is based on orthogonal polynomials.)

The rate at which an iterative method converges depends greatly on the spectrum of the coefficient matrix. Hence, iterative methods usually involve a second matrix that transforms the coefficient matrix into one with a more favorable spectrum. The transformation matrix is called a preconditioner. A good preconditioner improves the convergence of the iterative method, sufficiently to overcome the extra cost of constructing and applying the preconditioner. Indeed, without a preconditioner the iterative method may even fail to converge.

3.1 The why and how

The convergence rate of iterative methods depends on spectral properties of the coefficient matrix. Hence one may attempt to transform the linear system into one that is equivalent in the sense that it has the same solution, but that has more favorable spectral properties. A preconditioner is a matrix that effects such a transformation.

For instance, if a matrix M approximates the coefficient matrix A in some way, the transformed system M−1 Ax= M−1 b has the same solution as the original system Ax=b , but the spectral properties of its coefficient matrix M−1 A may be more favorable.

In devising a preconditioner, we are faced with a choice between finding a matrix M that approximates A, and for which solving a system is easier than solving one with A, or finding a matrix M that approximates A−1 , so that only multiplication by M is needed. The majority of preconditioners falls in the first category; a notable example of the second category will be discussed in §3.5.

4.1 Complex Systems

Conjugate gradient methods for real symmetric systems can be applied to complex Hermitian systems in a straightforward manner. For non-Hermitian complex systems we distinguish two cases. In general, for any coefficient matrix a CGNE method is possible, that is, a conjugate gradients method on the normal equations AH Ax= AH b , or one can split the system into real and complex parts and use a method such as GMRES on the resulting real nonsymmetric system. However, in certain practical situations the complex system is non-Hermitian but symmetric.

Complex symmetric systems can be solved by a classical conjugate gradient or Lanczos method, that is, with short recurrences, if the complex inner product (x,y) = x ¯ T y is replaced by (x,y) = xT y . Like the BiConjugate Gradient method, this method is susceptible to breakdown, that is, it can happen that xT x=0 for x≠0 . A look-ahead strategy can remedy this in most cases (see Freund [100] and Van der Vorst and Melissen [208]).

5.1 The Lanczos Connection

As discussed by Paige and Saunders in [168] and by Golub and Van Loan in [109], it is straightforward to derive the conjugate gradient method for solving symmetric positive definite linear systems from the Lanczos algorithm for solving symmetric eigensystems and vice versa. As an example, let us consider how one can derive the Lanczos process for symmetric eigensystems from the (unpreconditioned) conjugate gradient method.

A large body of numerical software is freely available 24 hours a day via an electronic service called Netlib. In addition to the template material, there are dozens of other libraries, technical reports on various parallel computers and software, test data, facilities to automatically translate FORTRAN programs to C, bibliographies, names and addresses of scientists and mathematicians, and so on. One can communicate with Netlib in one of a number of ways: by email, through anonymous ftp (netlib2.cs.utk.edu) or (much more easily) via the World Wide Web through some web browser like Netscape or Mosaic. The url for the Templates work is: http://www.netlib.org/templates/. The html version of this book can be found in: http://www.netlib.org/templates/Templates.html.

The BLAS give us a standardized set of basic codes for performing operations on vectors and matrices. BLAS take advantage of the Fortran storage structure and the structure of the mathematical system wherever possible. Additionally, many computers have the BLAS library optimized to their system. Here we use five routines:

1. SCOPY: copies a vector onto another vector

2. SAXPY: adds vector x (multiplied by a scalar) to vector y

3. SGEMV: general matrix vector product

4. STRMV: matrix vector product when the matrix is triangular

5. STRSV: solves Tx=b for triangular matrix T

Appendix C: Glossary

The back matter includes bibliography and index.