Incomplete LU Preconditioners for Conjugate-Gradient-Type Iterative Methods

Abstract

Summary In recent years, such conjugate-gradient-type methods as orthomin have been used very successfully with various preconditioners to solve the unsymmetric linear systems that arise in reservoir simulation. Here these successful iterative methods are combined with a new set of preconditioners that have been derived from some powerful algorithms in direct sparse Gaussian elimination. The SPARSPAK implementations of minimum degree and nested dissection algorithms have been modified to yield incomplete LU factorizations. Their applications to reservoir simulation problems is investigated here. Introduction This paper describes a general sparse-linear-equation solver called ILUPACK,1 which consists of three phases: reordering, incomplete factorization, and iterative solution. For each phase, a set of different algorithms is provided. Any possible combination of algorithms for the different phases can be used to obtain a solution method. Thus the ILUPACK linear equation solver provides a good experimental tool to explore the range of applications and the efficiency of a particular solution technique. On the basis of experimental results with ILUPACK, a special-purpose linear equation solver can be developed, tailored to the particular application at hand. ILUPACK is available on the MAINSTREAM EKS/VSP (Cray X-MP) service of Boeing Computer Services. There are no current plans to distribute the software because of its experimental nature. A new package with capabilities similar to ILUPACK is under development.2 This package will be tailored more closely to the architecture of current vector computers. The user interface of ILUPACK is based on a set of subroutines developed by George and Liu3 for SPARSAK. This user interface allows easy entry of general sparse matrices. No a priori restrictions are made with respect to a block or band structure of the matrix, or with respect to the number of unknowns per gridpoint. Thus the package can easily handle irregular grids, such as those arising from reservoirs with fault links, in reservoirs with refined grids around wells, or in reservoirs with coarser two-dimensional grids modeling a surrounding aquifer. ILUPACK provides the user with the capability of quickly evaluating new approaches to model difficult problems. It relieves the user of the requirement to develop special-purpose factorization modules for each such problem. On the other hand, the complete generality of ILUPACK has its price, because special properties of the underlying physical problem cannot be recovered from the purely algebraic treatment of the matrix problem. We believe, however, that the capability of solving very general problems outweighs this disadvantage. Reordering Methods The solution of sparse linear systems of equations of the formEquation 1 by direct methods has been an area of intensive research during the past 15 years. Most of the effort has been directed toward a combination of Gaussian elimination with some reordering of the equations and unknowns in Eq. 1. The goal is to obtain a permuted system,Equation 2 that can be solved more easily. Here A=PAQ, x=QTx, and b=Pb. Here P and Q are permutation matrices chosen, for example, such that the factorization of A incurs less fill-in than the factorization of A. For general nonsymmetric A, the permutations P and Q are determined during the course of the numerical factorization and usually depend on the numerical value of the entries in A. The situation becomes much easier when a is symmetric and positive definite. In this case, reorderins can be found on the basis of the zero/nonzero structure of A alone. The numerical values of the entries of A are irrelevant, because if A is positive definite, then so is A=PAPT, and a Choleski factorization, A=LLT, can always be computed. Because of this decoupling of ordering and factorization phase, the symmetric positive-definite case is easier than the general case, and the state of sparse direct methods for symmetric positive-definite problems is very satisfactory. The work by George and Liu4 can be regarded as the culmination of research in that area. Their sparse-matrix package, SPARSPAK, provides high-quality implementations of several algorithms. Of particular interest to us are some of their reordering algorithms. The idea of combining a general matrix reordering algorithm with an incomplete factorization is fairly new. This is the main topic of this report. A heuristic argument can be used to justify why a reordering algorithm should add to the benefits of incomplete factorizations. The goal of any reordering is to produce a permutation so that the permuted system (Eq. 2) has less fill-in than the original system (Eq. 1). The number of nonzeros in the factors of A is less than the number of nonzeros in the factors of A. Suppose an incomplete factorization of both A and A is performed and all fill-in is discarded - i.e., only an incomplete factorization based on the structure of A or A is computed. Then, because in the system of Eq. 2 a full factorization would have incurred less overall fill-in than in Eq. 1, the incomplete factor captures a larger share of information of the complete factor in Eq. 2 than in Eq. 1. In an extreme case, one could think of an example in which a reordering would result in no fill-in whatsoever. Then the incomplete factors of Eq. 2 would be identical to the complete factors of Eq. 2, whereas the incomplete factors of Eq. 1 still would be quite different from the complete factors of Eq. 1. Thus, for a comparatively small effort, a reordering may improve the quality of an incomplete factorization.