**Preface to the existing class notes
**

At the risk of mixing notation a little I want to discuss the general form of

iterative methods at a general level. We are trying to solve a linear system

situation where cost of direct solution (e.g. Gauss elimination) for matrix

One approach starts by writing

that a linear system in the form

rearrange the original linear equation to the form:

Introduce a series of approximations to the solution x starting with an initial guess

and proceeding to the latest available approximation . The next approximation is

obtained from the equation:

or

Understand that we never actually generate the inverse of E, just solve the equation for

To understand convergence properties of this iteration introduce the residual

. Using this definition and the definition of the iteration, we get the useful

expression:

or

If you know something about linear algebra this tells you that the absolute values of

eigenvalues of

Jacobi, Gauss-Seidel, SOR, and ADI methods.

Another approach is to choose a matrix

close to the identity matrix and the solution of an equation like

without too much effort. The revised problem can be cast either in the form

or

This use of M is referred to as matrix preconditioning. For any given iteration the error

can be approximated as

giving a new iteration form of

If we start with a guess of , then the iteration yields approximations which

are combinations of the vectors {b,

of vectors is called a Krylov subspace. Krylov solution methods apply the above

iteration philosophy, and modify the coefficients to the vectors in the Krylov subspace to

somehow minimize the residual or error associated with each approximation. This class

of iteration techniques includes the Conjugate Gradient (CG), MINRES, GMRES, QMR,

BiCGSTAB, and a host of related methods. They generally outperform the first class of

iteration by a significant margin.

Fluid flow and heat transfer codes generate sparse linear systems in two contexts.

The first is a nearly block tridiagonal system resulting from the equations in 1-D models.

Direct solution methods have worked well for these systems. However, the details of this

approach need reconsideration, if codes are adapted to massively parallel computers. The

second, and more problematic class of sparse matrices are generated by equations

modeling two- and three-dimensional regions. These systems are often only marginally

diagonally dominant, and hence pose a significant challenge to iterative solution

methods. Twenty years ago finding a method that dealt with these equations well was

extremely difficult. Now the biggest problem is sorting through a wide selection of

methods for the problem to find the most appropriate approach.

About fifteen years ago enough experience had been gained with matrix

preconditoning that variants on the Conjugate Gradient method had come into wide

use. Evolution of this methodology has continued with the introduction of several

variations on the basic algorithm. The most popular of these is currently Sonneveld's

conjugate gradient squared (CGS) algorithm . This class of methods has a rate of

convergence that is generally very good, but is not monotonic. Plots of residual versus

iteration count can show oscillations. A stabilized CGS method (CGSTAB) has been

introduced to mitigate the oscillatory behavior of residuals, but it does not guarantee

monotonically decreasing residuals.

Conjugate Gradient methods currently are losing favor to more general Krylov

subspace methods based largely on the GMRES algorithm . As with conjugate

gradient, this method is based on the construction of a set of basis vectors, and formally

will converge to the exact solution. Rapid convergence in the initial iterations requires

preconditioning of the matrix in both approaches. GMRES type algorithms have the

advantage that residuals decrease monotonically, and that the algorithms are generally

more robust. They have the disadvantage that they must store an additional basis vector

in the Krylov subspace for each iteration. The partial solution to this problem has been to

restart the solution algorithm after some number of iterations.

Providing a recommendation for a "best" solution algorithm is not currently

possible. In fact variability of algorithm performance with machine architecture and

problem type suggests that a "best" algorithm exists only in an average sense. Some

guidance can be obtained from recent publications. Soria and Ruel provide a

summary of the algorithms mentioned above, and comparisons of GMRES and CGS

using ILU or diagonal preconditioning to results of Broyden , Gauss-Seidel, and

Jacobi iterations. The authors clearly illustrate the value of good preconditioning, and

conclude that more than one solution procedure should be available in a CFD program.

However, they offer no advice on criteria for method selection in such a CFD program.

Chin and Forsyth provide more detailed information for
judging relative

performance of GMRES and CGSTAB with ILU preconditioning. For their problems,

CGSTAB was generally faster than GMRES, but they noted GMRES was not afflicted

with the occasional divide by zero observed in the CGSTAB algorithm. They also
note

(but don't document) the sensitivity of conclusions on relative performance of
these

methods to the quality of the initial guess for the problem solution.

The interaction between algorithm and machine architecture is always a popular

topic. Recent discussions and references can be found in articles by van Gijzen

(impact of vectorization on GMRES), Sturler and van der Vorst (GMRES and CG
on

distributed memory parallel machines), and Xu et al. (GMRES for parallel

machines). However, for production codes with a wide user base it is more
important to

first select a method that performs well without special vector or parallel
considerations.

This helps to minimize non-uniform behavior across platforms. Special vectorized
or

parallelized packages should be options that can be activated and checked after
initial

code validation on a new installation.

I would recommend several key considerations in choosing a linear system solver.

The starting point, of course, is a set of problems (probably only segments of
transients)

that are believed to be representative of the code's workload. Enough solution
packages

are publicly available that special programming of algorithms should not be
attempted for

initial tests. LASPack is one resource for testing combinations of preconditioners

and solution algorithms. Once a good (based on robustness and speed) iterative
solver

has been selected, the break-even point (in system size) should be determined
between

use of the iterative solution and an efficient direct sparse matrix solver.
Frequently both

iterative and direct solution packages should be included in a simulation code
with an

internal check on the geometry of the mesh to choose between the methods.

If you are willing to wade through lots of Lemmas and Theorems, I recommend

that you read a recent (1997) work by Anne Greenbaum , summarizing the state
of

the art at that time. She includes a vary extensive list of references to the
literature on

iterative methods.