Ordering strategies and related techniques to overcome the trade-off between parallelism and convergence in incomplete factorizations

Abstract

This paper is concerned with the parallel implementation of the incomplete factorization preconditioned iterative method. Although the use of such parallel ordering as multicolor ordering may increase parallelism in factorization, it often slows convergence when used in the preconditioned method, and thus may offset the gain in speed obtained with parallelization. Further, the higher the parallelism of an ordering, the slower the convergence; the lower the parallelism, the faster the convergence. This well-known trade-off between parallelism and convergence is well explained by the property of compatibility, the level of which can be clearly seen when ordering is presented in graph form (S. Doi, A. Lichnewsky, A graph-theory approach for analyzing the effects of ordering on ILU preconditioning, INRIA report 1452, 1991). In any given method, the fewer the incompatible local graphs in an ordering (i.e., the lower the parallelism), the faster the convergence (S. Doi, Appl. Numer. Math. 7 (1991) 417–436; S. Doi, in: T. Nodera (Ed.), Advances in Numerical Methods for Large Sparse Sets of Linear Systems, 7 Keio University, 1991). An ordering with no incompatible local graphs, for example, such as that implemented on vector multiprocessors by using the nested dissection technique, will have excellent convergence, but its parallelism will be limited (S. Doi, A. Lichnewsky, Int. J. High Speed Comput. 2 (1990) 143–179). To attain a better balance, a certain degree of incompatibility is necessary. In this regard, increasing the number of colors in multicolor ordering can be a useful approach (S. Fujino, S. Doi, in: R. Beauwens (Ed.), Proceeding of the IMACS Internation Symposium on Iterative Methods in Linear Algebra, March 1991; S. Doi, A. Hoshi, Int. J. Comput. Meth. 44 (1992) 143–152). Two related techniques also presented here are the overlapped multicolor ordering (T. Washio, K. Hayami, SIAM J. Sci. Comput. 16 (1995) 631–650), and a fill-in strategy selectively applied to incompatible local graphs. Experiments conducted with an SX-5/16A vector parallel supercomputer show the relative effectiveness of increasing the number of colors and also of using this approach in combination with overlapping and with fill-ins.

Introduction

A class of the most effective methods for solving a system of large sparse linear equations is the preconditioned iterative method, which can be defined as a basic iterative method applied to a preconditioned system M⁻¹Au=M⁻¹b (where M is the preconditioner) instead of being applied to the original system Au=b. The objective of the preconditioning is to reduce the condition number (or to cluster the eigenvalues) of the original system so as to reach an approximate solution with fewer iterations. This requires solving the system Mv=g. Hence, a good preconditioner needs to satisfy two requirements: it should be a good enough approximation of A to produce a reduced condition number (or better clusters of eigenvalues), and the system Mv=g should be much easier to solve than the original system. Further, increases in the parallelism of recent computer architectures have led to a new requirement: the solution of Mv=g must also have high enough parallelism to be mapped naturally on the computer to be used.

One type of commonly used preconditioning is represented by the incomplete factorization M=LU, where L and U are lower and upper triangular matrices of a sparsity structure similar to A, and are produced by neglecting certain fill-ins in the Gaussian elimination [19]. One way to attain parallelism in the solution of Mv=LUv=g (in practice, this solution is obtained first by solving Lw=g (forward substitution) and then by solving Uv=w (backward substitution)) is to reorder the system Au=b and to reconstruct an incomplete factorization for the reordered system. Actually, a number of techniques developed in the past do this [1], [3], [5], [7], [10], [12], [13], [14], [17], [18], [21]. In this regard, one possibility for this approach is to use the well-known red–black ordering. Solving Lw=g, where w has a red–black ordering can be performed simultaneously for half the unknowns in v; this is also true of the solving of Uv=w. In this case, however, a fundamental difficulty is presented by the trade-off problem between parallelism and convergence: in a system to which the incomplete factorization preconditioned iterative method applied, higher parallelism in the ordering will result in slower convergence, and higher convergence will require lower parallelism. This trade-off problem was first reported by Duff and Meurant [10] on the basis of intensive numerical tests using many orderings [10]. Since then, many attempts have been made to answer the following two fundamental questions: why there is a trade-off between parallelism and convergence in incomplete factorizations and how can the trade-off be overcome [2], [3], [6], [10], [11], [16], [20]? Before we give our own answers to these questions, it is worth repeating the remarks of Duff and Meurant [10] which are still valuable:

• (DM1) The rate of convergence is almost directly related to the norm of the remainder matrix R=M−A, but not related to the number of fill-ins dropped.

• (DM2) It appears that the harder the problem (discontinuous coefficients, anisotropy, etc.) is, the more important the ordering for the incomplete factorization is.

• (DM3) A single level of fill-in changes the relative performance of the different ordering schemes; for example, mind, rb, and altd (parallel orderings) do somewhat better when some fill-in is allowed in L. One reason for this is that many of the first-level fill-ins for these orderings are quite large, unlike row, for example, where the fill-ins rapidly decrease in value.

Kuo and Chan [16] proved that the condition number of a preconditioned system based on a red–black ordering is only about 1/4 of that of an unpreconditioned system, showing no asymptotic improvement no matter how small grid size h becomes [16]. Eijkhout [11] used an analysis of the infinity norm of the remainder matrix R=M−A to create a criterion by which orderings for no-fill factorization might be placed in one of two categories: whether or not an ordering contains nodes that are eliminated before both neighbors in one direction are eliminated [11].

Doi and Lichnewsky [6] attempted to answer the first question “why is there a trade-off(?)” by utilizing graph representations of orderings [6]. They showed that when an ordering graph contains incompatible (local) graph, the components of the remainder matrix R corresponding to the incompatible graph are significantly larger than the components corresponding to any single compatible graph. Furthermore, while components corresponding to a compatible graph converge to zero with dynamic increase or decrease in the value of model equation parameters, (e.g., grid aspect ratio, anisotropy, etc.) components corresponding to an incompatible graph do not. Another interesting observation is that the components of r^k=b−Au^k (the residual vector after k iterations with a preconditioned iterative method) corresponding to incompatible local graphs have significantly larger values than those corresponding to compatible local graphs.¹ An ordering which has no incompatible graphs is called compatible. Doi and Lichniewsky also confirmed that the trade-off observed by Duff and Meurant could be explained very well by means of the quality of compatibility. With respect to parallel computing, the incompatible local graph corresponds to a potential starting point for forward substitution by L. This may explain the trade-off; that is, an ordering having higher parallelism has more incompatible local graphs that break the compatibility of the ordering and thus decrease the rate of convergence. It is also interesting to see that compatibility is equivalent to Eijkhout's criterion, although his study was performed independently. Doi and Lichnewsky [5] implemented several incomplete factorizations based on compatible orderings on a Cray-2 shared memory vector multiprocessor using the nested dissection technique, and showed that these orderings can be effectively mapped onto vector multiprocessors [5]. Doi [3], [4] did a more quantitative study of this trade-off problem by using two parameterized orderings, parallel block PB(m) and parallel diagonal PD(m) (where m is a parameter that controls the degree of parallelism), and showed that the rate of convergence degrades rather smoothly as the degree of parallelism increases [3]. A pragmatic solution designed to strike a better balance between parallelism and convergence is the use of 50 or more colors in multicolor ordering, rather than the conventionally used four or eight colors [4], [7], [13]. Such a large-numbered multicolor ordering technique attains high performance on an SX-3/14 vector supercomputer with minor degradation in convergence [7].

Some positive techniques have also been proposed to compensate for the significantly larger components of R that corresponding to the incompatible graphs which occasionally appear in large-numbered multicolor incomplete factorization. Washio and Hayami [24], for example, introduced the overlapping technique. The key idea is to compensate for those components by repeating parts of the forward and backward substitutions [24]. It is also possible to apply the fill-in technique to incompatible graphs. Since those appear only very rarely in large numbered multicolor ordering, this fill-in does not significantly reduce parallelism, and the cost of the fill-ins is sufficiently low.

Section 2 of this paper introduces some fundamental definitions, including equivalence and compatibility in ordering. Section 3 presents analysis that answers the question “why is there a trade-off?” Section 4 introduces ordering strategies and related techniques; those are in some sense answers to the question “how can we overcome the tradeoff?” Section 5 presents the results of numerical experiments conducted on an SX-5 vector parallel supercomputer. Our results show that a large-numbered multicolor MILU–BiCGSTAB [23] method, combined with either overlapping or selective fill-ins, can solve some discrete convection–diffusion equations, discretized on a 199³ grid, within a time on the order of 1 s, where sustained speed is about 40 Gflops on an SX-5/16A (16 CPUs, peak speed of 128 Gflops).