Theoretical Div., Mail stop B-214 Los Alamos National Laboratory Los Alamos, N.M. 87545

Abstract

We present a novel linear solver which works well for large systems obtained from discretizing a PDE. In particular, it is robust and the computational effort scales as N log(N) where N is the number of equations. The algorithm is based on a wavelength decomposition which 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 novel choice of the preconditioning matrix. The preconditioner is an approximate inverse of the linear operator. It is constructed from the inverse of the coarse grained linear operator and from smoothing operators which 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 its inverse is not known explicitly, the conjugate gradient 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 is terminated 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 extra 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. For simple problems, such as arise from the Poisson equation, later iterations on the finer scales need fewer iterations on the coarser scales and the computational effort is proportional to N. For more complicated problems, operations performed on the finest scale still dominate but the computational effort scales as N log(N). We have tested our solver on the porous flow equation. On a workstation we have solved problems on grids as large as 1000x1000. The algorithm works well even when the permeability tensor has spatial variations exceeding a factor of 10^6.