Algorithms for Solving Linear Systems

In this paper we are interested in methods for solving large sparse linear systems that arise from the discretization of partial differential equations, denoted by

$$\begin{eqnarray*}Ax=b, \end{eqnarray*}$$

where $A$ is an $N\times N$ matrix. In this section we give a brief overview of the methods we have used in our scaling studies.

One of the leading families of methods for linear systems are iterative solvers known as Krylov subspace methods [6], [8], [9], [12] . These methods are characterized by the subspaces in which the iterates lie. The $j$-th Krylov subspace for a given matrix $A$ and vector $r_0$ is

$$\begin{eqnarray*}K_j(A,r_0) = \span\{r_0, Ar_0, A^2 r_0,\dots, A^{j-1}r_0\}. \end{eqnarray*}$$

The $j$-th iterate of a Krylov subspace method lies in the shifted affine subspace

$$\begin{eqnarray*}x_j \in x_0 + K_j(A,r_0). \end{eqnarray*}$$

Different Krylov subspace methods differ in the criteria they use in selecting a vector in the subspace. The Conjugate Gradient method picks $x_j$ to be the vector that minimizes the $A$-norm of the error over the Krylov subspace. This method can only be used for problems where the matrix is symmetric and positive definite.

A number of Krylov subspace methods have been developed for the case when the matrix is not symmetric, positive definite. The GMRES method [11] picks $x_j$ to be the vector that minimizes the two norm of the residual over the Krylov subspace. In its original form GMRES requires one additional vector of storage for each iteration, for the Arnoldi process, that applies Gram-Schmidt orthogonalization to the Krylov basis. In practice, GMRES is usually implemented in its restarted form, GMRES(k), with a fixed number of vectors stored. Another popular Krylov subspace method for non-symmetric problems is the BiCGStab (stabilized BiConjugate Gradient) method [13].

The main operations in a Krylov subspace method are (i) matrix-vector products (ii) dot products and norm computations, and (iii) vector updates.

The convergence rate of a Krylov subspace method can be improved by using the method in conjunction with a preconditioner, which we denote by $M$. The iterative method is applied to the preconditioned system

$$\begin{eqnarray*} M^{-1} A x = M^{-1} b. \end{eqnarray*}$$

Applying the preconditioner to a vector $r$ is symbolically denoted by $z = M^{-1} r$. In selecting a preconditioner one needs to balance the cost of applying the preconditioner and the improvement in convergence obtained by preconditioning (which is related to how well $M$ approximates $A$). Some of the simplest preconditioners are Jacobi and block Jacobi. Other preconditioners are ILU (incomplete LU factorization) and approaches based on domain decomposition, such as Additive Schwartz.

Another method that can be used on its own, or as a preconditioner in a Krylov subspace method, is the multigrid method [7]. In its usual (geometric) form, multigrid works with a discretization of the PDE on a hierarchy of progressively coarser grids. The idea is to apply a few iterations of a relaxation method at a given level to get an approximation whose error is smooth and can be represented on a coarser grid. Then the residual of this approximation is restricted to the next coarser grid, and the process repeated until the coarsest grid is reached. Here a coarse grid solve is performed, and the solution (the coarse grid correction) is interpolated up to the finer grid. After adding the coarse grid solution to the existing approximation, a few more iterations of the smoother are applied.

The user of a multigrid code needs to make a number of choices to get a specific form of the method. These choices include the smoother, the number of pre- and post-smoothing iterations, the restriction and interpolation operators, coarse grid matrices for each grid, the coarse grid solver, and the cycling scheme. The preconditioned Conjugate gradient method requires the preconditioning matrix to be symmetric. Since the multigrid iteration matrix is not symmetric (unless special choices are made to ensure symmetry), we use non-symmetric Krylov subspace methods in conjunction with multigrid preconditioners.