Computing Systems in Engineering Vol. 1, Nos 2-4, pp. 183-195, 1990 0956-0521/90 $3.00+0.00 Printed in Great Britain. © 1990 Pergamon Press plc A PARALLEL HYBRID SPARSE LINEAR SYSTEM SOLVER K. GALL|VAN,t A. SAMEHt and Z. ZLATEV~ tCenter for Supercomputing Research and Development, University of Illinois at Urbana-Champaign, 305 Talbot Laboratory, 104 South Wright Street, Urbana, IL 61801, U.S.A. ~tAir Pollution Laboratory, Danish Agency of Environmental Protection, Risoe National Laboratory, DK-4000 Riskilde, Denmark (Received 27 April 1990) Abstraet-42onsider the system Ax = b, where A is a large sparse nonsymmetric matrix. It is assumed that A has no sparsity structure that may be exploited in the solution process, its spectrum may lie on both sides of the imaginary axis and its symmetric part may be indefinite. For such systems direct methods may be both time consuming and storage demanding, while iterative methods may not converge. In this paper, a hybrid method, which attempts to avoid these drawbacks, is proposed. An L U factorization of A that depends on a strategy that drops small non-zero elements during the Gaussian elimination process is used as a preconditioner for conjugate gradient-like schemes, ORTHOMIN, GMRES and CGS. Robustness is achieved by altering the drop tolerance and recomputing the preconditioner in the event that the factorization or the iterative method fails. If after a prescribed number of trials the iterative method is still not eonvergent, then a switch is made to a direct solver. Numerical examples, using matrices from the Harwell-Boeing test matrices, show that this hybrid scheme is often less time consuming and storage de- manding; than direct solvers, and more robust than iterative methods that depend on preconditioners that depend .an classical positional dropping strategies. I, THE HYBRID ALGORITHM Consider the system of linear algebraic equations Ax = b, where A is a nonsingular, large, sparse and nonsymmetric matrix. We assume also that matrix A is generally sparse (i.e. it has neither any special property, such as symmetry and/or positive definite- ness, nor any special pattern, such as bandedness, that can be exploited in the solution of the system). Solving such linear systems may be a rather difficult task. This is so because commonly used direct methods (sparse Gaussian elimination) are too time consuming, and iterative methods whose success depends on the matrix having a definite symmetric part or depends on the spectrum lying on one side of the imaginary axis are not robust enough. Direct methods have the advantage that they normally produce a sufficiently accurate solution, although a direct estimation of the accuracy actually achieved requires additional work. On the other hand, when iterative methods converge sufficiently fast, they require computing time that is several orders of magnitude smaller than that of any direct method. This brief comparison of the main properties of direct methods and iterative methods for the problem at hand shows that the methods of both groups have some advantages and some dis- advantages. Ttlerefore it seems worthwhile to design methods that combine the advantages of both groups, while minimizing their disadvantages. Throughout we assume that sparse Gaussian elimination is the direct method chosen for the solution of Ax = b. It is well known that this is the best choice in the case where A is large and generally sparse; see for example Ref. 1 or 2. The arithmetic operations during stage k(k = 1,2...n- 1 ) of Gaussian elimination are carried out by the formula a!k+~ ~(k}__.(k)t.(Ap~--l.(k) l) q = ¢til ~ik ~,Ukk ] tSkj , where i=k+l, k+2...n, j=k+l, k+2...n, while (k) aii = a,~ are the elements of matrix A. It is clear that = {,~ a},} ~ vanish, then a if a}~) 0 while neither ai, nor new non-zero element, fill-in, is created in position (i, j). Unfortunately, fill-in does occur when large sparse matrices are factored by Gaussian elimination and this method becomes rather expensive (in terms of time and storage requirements) when many fill-ins are introduced. Therefore, reducing the number of fill-ins is one of the main tasks during the develop- ment of sparse matrix cocles, at least on sequential and vector computers. Such minimization of fill-in is achieved by adopting a suitable pivoting strategy; see for example Refs 3-5. The amount of fill-in may be large, however, even when a good pivoting strategy is adopted. For such systems direct methods may lose competitiveness with iterative methods if a con- vergent method can be found for the system. The successful use of iterative methods often depends upon the effective use of preconditioning. Precon- ditioning a system Ax = b involves using the iterative method to solve the related system M -1 Ax = M-~b, where the preconditioned M is easily invertible and M -1 A ~ I. The choice of M is an art in itself and de- pends on the iterative method used, the application from which the system arises, and, for high-per- formance machines, the architecture on which the algorithm is to execute. 6 For the purpose of obtaining an approximate faetorization of A for preconditioning, while main- 183