Multi-Scale Linear Solvers for Very Large Systems Derived from PDEs

        Klaus Lackner (ksl@lanl.gov)
                Theoretical Division, Mail Stop B-216,
                Los Alamos National Laboratory,
                Los Alamos, NM 87544
        
        Ralph Menikoff (rtm@lanl.gov)
                Theoretical Division, Mail Stop B-214,
                Los Alamos National Laboratory,
                Los Alamos, NM 87544
                 
ABSTRACT

We present a novel linear solver that works well for large systems obtained
from discretizing PDEs.  It is robust and, for the examples we studied, the
computational effort scales linearly with the number of equations.  The
algorithm is based on a wavelength decomposition that combines conjugate
gradient, multi-scaling and iterative splitting methods into a single
approach.  On the surface, the algorithm is a simple preconditioned conjugate
gradient with all the sophistication of the algorithm in the choice of the
preconditioning matrix.  The preconditioner is a very good approximate inverse
of the linear operator.  It is constructed from the inverse of the coarse
grained linear operator and from smoothing operators that are based on an
operator splitting on the fine grid.  The coarse graining captures the long
wavelength behavior of the inverse operator while the smoothing operator
captures the short wavelength behavior.  The conjugate gradient iteration
accounts for the coupling between long and short wavelengths.  The coarse
grained operator corresponds to a lower resolution approximation to the PDEs.
While the coarse grained inverse is not known explicitly, the algorithm only
requires that the preconditioner can be a applied to a vector.  The coarse
inverse applied to a vector can be obtained as the solution of another
preconditioned conjugate gradient solver that applies the same algorithm to
the smaller problem.  Thus, the method is naturally recursive.  The recursion
ends when the matrix is sufficiently small for a solution to be obtained
efficiently with a standard solver.  The local feedback provided by the
conjugate gradient step at every level makes the algorithm very robust.  In
spite of the effort required for the coarse inverse, the algorithm is
efficient because the increased quality of the approximate inverse greatly
reduces the number of times the preconditioner needs to be evaluated.  A
feature of the algorithm is that the transition between coarse grids is
determined dynamically by the accuracy requirement of the conjugate gradient
solver at each level.  Typically, later iterations on the finer scales need
fewer iterations on the coarser scales and the computational effort is
proportional to N rather than NlogN, where N is the number of equations.  We
have tested our solver on the porous flow equation.  On a workstation we have
solved problems on grids ranging in dimension over 3 orders of magnitude, from
10<sup>3</sup> to 10<sup>6</sup>, and found that the linear scaling holds.
The algorithm works well, even when the permeability tensor has spatial
variations exceeding a factor of 10<sup>9</sup>.