A Multilevel Dual Reordering Strategy for Robust Incomplete LU Factorization of Indefinite Matrices.
A dual reordering strategy based on both threshold and graph reorderings is introduced to construct robust incomplete LU (ILU) factorization of indefinite matrices. The ILU matrix is constructed as a preconditioner for the original matrix to be used in a preconditioned iterative scheme. The matrix is first divided into two parts according to a threshold parameter to control diagonal dominance. The first part with large diagonal dominance is reordered using a graph-based strategy, followed by an ILU factorization. A partial ILU factorization is applied to the second part to yield an approximate Schur complement matrix. The whole process is repeated on the Schur complement matrix and continues for a few times to yield a multilevel ILU factorization. Analyses are conducted to show how the Schur complement approach removes small diagonal elements of indefinite matrices and how the stability of the LU factor affects the quality of the preconditioner. Numerical results are used to compare the new preconditioning strategy with two popular ILU preconditioning techniques and a multilevel block ILU threshold preconditioner.
Introduction
This paper is concerned with reordering strategies used in developing robust preconditioners based
on incomplete LU (ILU) factorization of the coefficient matrix of sparse linear system of the form
where A is an unstructured matrix of order n. In particular, we are interested in ILU preconditioning
techniques for which A is an indefinite matrix; i.e., a matrix with an indefinite symmetric part.
Indefinite matrices arise frequently from finite element discretizations of coupled partial differential
equations in computational fluid dynamics and from other applications.
Technical Report 285-99, Department of Computer Science, University of Kentucky, Lexington, KY, 1999. This
work was supported in part by the University of Kentucky Center for Computational Sciences and in part by the
University of Kentucky College of Engineering.
y E-mail: jzhang@cs.uky.edu. URL: http://www.cs.uky.edu/~jzhang.
ILU preconditioning techniques have been successful for solving many nonsymmetric and indefinite
matrices, despite the fact that their existence in these applications is not guaranteed. However,
their failure rates are still too high for them to be used as blackbox library software for solving general
sparse matrices of practical interests [9]. In fact, the lack of robustness of preconditioned iterative
methods is currently the major impediment for them to gain acceptance in industrial applications,
in spite of their intrinsic advantage for large scale problems.
For indefinite matrices, there are at least two reasons that make ILU factorization approaches
problematic [9]. The first problem is due to small or zero pivots [23]. Pivots in an indefinite matrix can
be arbitrarily small. This may lead to unstable and inaccurate factorizations. In such cases, the size
of the elements in the LU factors may be very large and these large size elements lead to inaccurate
factorization. The second problem is due to unstable triangular solves [18]. The incomplete factors
of an indefinite matrix are usually not diagonally dominant. An indication of unstable triangular
solves is when kL are extremely large while the offdiagonal elements of L and U are
reasonably bounded. Such problems are usually caused by very small pivots. They may sometimes
happen without a small pivot. A statistic, condest, was introduced by Chow and Saad [9] to measure
the stability of the triangular solves. It is defined to be k(LU) \Gamma1 ek1 , where e is a vector of all ones.
This statistic is useful when its value is very large, e.g., in the order of 10 15 .
Small pivots are usually related small or zero diagonal elements. It can be argued that by
restricting the magnitude of the diagonal elements, we may be able to alleviate, if not eliminate,
these two problems of ILU factorizations to a certain degree. Such restrictions can be seen in the
form of full or partial pivoting strategies in Gaussian elimination. In ILU factorization, column
pivoting strategy has been implemented with Saad's ILUT, resulting in an ILUTP techniques [32].
However, ILUTP has not always been helpful in dealing with nonsymmetric matrices [3, 9]. As
Chow and Saad pointed [9], a poor pivoting sequence can occasionally trap a factorization into a zero
pivot, even if the factorization would have succeeded without pivoting. In addition, existing pivoting
strategies for incomplete factorization cannot guarantee that a nonzero pivot can always be found,
unlike the case with Gaussian elimination [9].
Another obvious strategy of dealing with small pivots is to replace them by a larger value. The
ILU factorization can continue and the resulting preconditioner may be well conditioned. In such a
way, the ILU factorization is said to be stabilized. However, this strategy alters the values of the
matrix and the resulting preconditioner may be inaccurate. Thus, the choice of the replacing value
for the small pivots is critical for a good performance and a good choice is usually problem dependent
[23]. Too large a value will result in a stable but inaccurate factorization; too small a value will result
in an unstable factorization. A similar strategy is to factor a shifted matrix A + ffI , where ff is a
positive scalar so that A+ ffI is well conditioned [27, 44]. Such a strategy too obviously has a tradeoff
between stable and accurate factorization. For more studies on the stability of ILU factorizations,
we refer to [19, 29, 42, 13, 45].
It is also possible to reorder the rows of the matrix so that their diagonal dominance in a certain
sense is in decreasing order. In this way, small pivots are in the last rows of the matrix and may not
be used in an ILU factorization. This strategy also has some problems since the values of the pivots
are modified in an unpredictable way, small pivots may still affect the ILU factorization. In addition,
the effect of standard reordering schemes applied to general nonsymmetric sparse matrices is still an
unsettled issue [17, 24, 43].
This paper follows the above idea of putting the rows with small diagonal elements to the last
few rows. However, these small diagonal elements will never be used in the ILU factorization. Instead,
these rows form the rows of a Schur complement matrix and the values of the diagonal elements are
modified in a systematic way. This process is continued for a few times until all small diagonal
elements are removed; or until the last Schur complement matrix is small enough that a complete
pivoting strategy can be implemented inexpensively. With this reordering strategy, we can expect to
obtain a stable and accurate ILU factorization. We also implement a graph based reordering strategy
(minimum degree algorithm) to reduce the fill-in amount during the stable ILU factorization.
This paper is organized as follows. The next section introduces a dual reordering strategy
based on both the values and the graph of the matrix. Section 3 discusses a partial ILU factorization
technique to construct the Schur complement matrix implicitly. Section 4 gives analyses on the values
of the diagonal elements of the Schur complement matrix and shows how the stability of the LU factor
affects the quality of a preconditioner. Section 5 outlines the multilevel dual reordering algorithm.
Section 6 contains numerical experiments. Concluding remarks are included in Section 7.
Reordering Strategy
Most reordering strategies are originally developed for the direct solution of sparse matrices based on
Gaussian elimination. They are mainly used to reduce the fill-in elements in the Gaussian elimination
process or to extract parallelism from LU factorizations [15, 22]. They have also been used in ILU
preconditioning techniques for almost the same reasons [16, 20, 30]. Various reordering strategies were
first studied for preconditioned conjugate gradient methods, i.e., for the cases where the matrix is
symmetric positive definite [1, 4, 5, 10, 11, 26, 31]. They were then extended for treating nonsymmetric
problems [2, 7, 12, 14]. Most of these strategies are based on the adjacency graph but not on the
values of the matrices. They are robust for general sparse matrices only if used with suitable pivoting
strategies, which are based on the values of the matrices, to prevent unstable factorizations. Hence,
reordering strategies based on matrix values are needed to yield robust stable ILU factorizations.
Such an observation has largely been overlooked in ILU techniques for some time, partly because the
early ILU techniques were mainly developed to solve sparse matrices arising from finite difference
discretizations of partial differential equations [28]. In such cases, the diagonal elements of the
matrices usually have nonzero values.
In this paper, we introduce a dual reordering strategy for robust ILU factorization for solving
general sparse indefinite matrices. To this end, we first introduce a strategy to determine the row
diagonal dominance of a matrix. 1 We actually compute a certain measure to determine the relative
strength of the diagonal element with respect to a certain norm of the row in question. Algorithm 2.1
is an example of computing a diagonal dominance measure for each row of the matrix and was
originally introduced in [41] as a diagonal threshold strategy in a multilevel ILU factorization.
1 The reference to row diagonal dominance is due to the assumption that our matrix is stored in a row oriented
format, such as in the compressed sparse row format [34]. The proposed strategy works equally well if the matrix is
stored in a column oriented format with the reference to column diagonal dominance.
Algorithm 2.1 Computing a measure for each row of a matrix.
1. For do
2. r
3. If r i 6= 0, then
4. ~ t
5. End if
6. End do
7.
8. For do
9.
10. End do
In Line 2 of the Algorithm 2.1 the set defined as
i.e., the nonzero row pattern for the row i. A row with a small absolute diagonal value will have a
small t i measure. A row with a zero diagonal value will have an exact zero t i measure.
denote the adjacency graph of the matrix A, where is the
set of vertices and K is the set of edges. Let (v denote an edge from vertex v j to vertex v k . Since
a node in the adjacency graph of a matrix corresponds to a row of the matrix, we will use the term
node and row of a matrix interchangeably. Given a diagonal threshold tolerance ffl ? 0, we divide the
nodes of A into two parts, V 1 and V 2 , such that
It is obvious that
For convenience, we assume that a symmetric permutation is performed so that the nodes in V 1
are listed first, followed by the nodes in V 2 . Since the nodes in V 1 are "good" for ILU factorization
in terms of stability, we may further improve the quality of the ILU factorization by implementing
a graph based reordering strategy. The following minimum degree reordering algorithm is just one
example of such graph based reordering strategies to reduce the fill-in elements in the ILU factorization

We denote by deg(v i ) the degree of the node v i , which equals the number of nonzero elements
of the ith row minus one; i.e.,
The set of the degrees of the rows of the matrix A can be conveniently computed when Algorithm 2.1
is run to compute the diagonal dominance measure of A. For example, in Line 2 of Algorithm 2.1,
the number of nonzero elements of the ith row will be counted.
After the first reordering based on the threshold tolerance ffl, we perform a second reordering
based on the degrees of the nodes. But the second reordering is only performed with respect to the
nodes in V 1 . To be more precise, we reorder the nodes in V 1 in a minimum degree fashion; i.e., the
nodes with smaller degrees are listed first, those with larger degrees are listed last. After the two
steps of reorderings, we have
are the permutation matrices corresponding to the threshold tolerance reordering
and the minimum degree reordering, respectively. We use P g here to emphasize that it is just a graph
based reordering strategy, and is not necessarily restricted to the minimum degree reordering. Other
graph based reordering strategies such as the Cuthill-McKee or reverse Cuthill-McKee algorithms [25]
may be used to replace the minimum degree strategy. But their meaning may be slightly changed
since not all neighboring nodes of a node in V 1 belong to V 1 , some of them may be in V 2 . For
simplicity, we use A to denote both the original and the permuted matrices in the sequel so that the
permutation matrices will no longer appear explicitly. We also refer to the two reordering strategies
as threshold reordering and graph reordering for short.
3 Partial ILU Factorization
An incomplete LU factorization process with a double dropping strategy (ILUT) is first applied to
the upper part (D F ) of the reordered matrix A in (2). The ILUT algorithm uses two parameters
and - to control the amount of fill-in elements caused by the Gaussian elimination process and is
described in detail in [32]. ILUT builds the preconditioning matrix row by row. For each row of
the LU factors, ILUT first drops all computed elements whose absolute values are smaller than -
times the average nonzero absolute values of the current row. After an (incomplete) row is computed,
ILUT performs a search with respect to the computed current row such that the largest p elements
in absolute values are kept, the rest nonzero elements are dropped again. Thus the resulting ILUT
factorization has at most p elements in each row of the L and U parts. The use of a double dropping
strategy ensures that the memory requirement be met. It is easy to see that the total storage cost
for ILUT is bounded by 2pn for a matrix of order n.
The ILUT process is continued to the second part of the matrix A in (2) with respect to the
C) submatrix. However, the elimination process is only performed with respect to the columns in
E, and linear combinations for columns in C are performed accordingly. In other words, the elements
corresponding to the C submatrix are not eliminated. Such a process is called a partial Gaussian
elimination or a partial LU factorization in [38]. Note that, due to the partial Gaussian elimination,
all rows in the (E C) submatrix can be processed independently (in parallel). This is because all
nodes in the E submatrix that are to be eliminated use only the computed (I)LU factorization of
the (D F ) part. Note also that the diagonal values of the rows of the C submatrix are never used
as pivots. It can be shown [38] that such a partial Gaussian elimination process modifies C into the
(incomplete) Schur complement of A. In exact arithmetic, C would be changed into
where LU is the standard LU factorization of the D submatrix. Hence, this method constructs the
Schur complement indirectly, in contrast to some alternative methods, e.g., the BILUM preconditioner
in [37], in which the Schur complement is constructed explicitly by matrix-matrix multiplications.
The partial ILU factorization process just described yields a block LU factorization of the matrix
A of the form ' D F
where I and 0 are generic identity and zero matrices, respectively. If the factorization is exact and if
we can solve the Schur complement matrix A 1 , the solution of the original linear system (1) can be
found by a backward substitution. This process is similar to the sparse direct solution method based
on one step cyclic reduction technique [22].
The partial ILU factorization process is the backbone of a domain based multilevel ILU pre-conditioning
technique (BILUTM) described in [38]. Such an ILU factorization with a suitable block
independent set ordering yields a preconditioner (BILUTM) that is highly robust and possesses high
degree of parallelism. However, in this paper, the parallelism due to block independent set ordering is
not our concern, we restrict our attention to the robustness of multilevel ILU factorization resulting
from removing small pivots.
We can heuristically argue that the ILU factorization resulting from applying the above partial
ILU factorization to the reordered matrix is likely to be more stable than that would be generated
by applying ILUT directly to the original matrix. This is because the factorization is essentially
performed with respect to the nodes in V 1 that have a relatively good diagonal dominance. The
partial ILU factorization with respect to the nodes in V 2 never needs to divide any pivot elements.
So there is no reason that large size elements should be produced.
As remarked previously, if we can solve the Schur complement matrix A 1 in (3) to a certain
degree, we can develop a two level preconditioner for the matrix A. An alternative is based on the
observation that A 1 is another sparse matrix and we can apply the same procedures to A 1 that have
been applied to A to yield an even smaller Schur complement A 2 . This is the philosophy of multilevel
ILU preconditioning techniques developed in [33, 37, 38]. However, for this moment, we only discuss
the possible construction of a two level preconditioner.
A two level preconditioner. The easiest way to construct a two level preconditioner is to apply
the ILUT factorization technique to the matrix A 1 . One question will be naturally asked: is the
ILUT factorization more stable when applied to A 1 than when applied to A?
Notice that since the nodes with good diagonal dominance have all been factored out, we tend
to think that the nodes of A 1 are not good for a stable ILUT factorization. This may not always be
true, since the measure of diagonal dominance computed in Algorithm 2.1 is relatively to a certain
norm of the row in question. We need to examine relative changes in size of the diagonal value when
a node is considered as a node in A and when it is considered as a node in A 1 .
4 Analyses
Diagonal submatrix D. For the easy of analysis, unless otherwise indicated explicitly, we assume
that the partial LU factorization described above is exact; i.e., no dropping strategy is enforced. We
also assume that, in the reordered matrix, the D submatrix is diagonal. Such a reordering can be
achieved by an independent set search as in a multielimination strategy of Saad [33, 39]. Thus, the
factorization (4) is reduced to
ED
'' D F
We now assume that all indices are local to individual submatrices. In other words, when we say the
ith row of the matrix F , we mean the ith row of the submatrix F , not the ith row of the matrix A,
original or permuted. For convenience we assume that D is of dimension m and A 1 is of dimension
m. We also use the notations:
It can be shown [22, 38] that, with the partial LU factorization without dropping, an arbitrary
element of the Schur complement matrix A 1 is:
Since we assume that the nodes with large diagonal dominance measure are in V 1 and the nodes in
have small or zero diagonal dominance measure, we are interested in knowing how the diagonal
value of a node of A may change when it becomes a node in A 1 .
row i
F
e ik
rowk
l
l

Figure

1: An illustration of the partial LU factorization to eliminate e ik in the E submatrix.
The following proposition is obvious from Equation (6) and from Figure 1.
Proposition 4.1 If either the jth column of the submatrix F or the ith row of the submatrix E is a
zero vector, then s
Definition 4.2 A node v i of the vertex set V is said to be independent to a subset V I of V if and
only if
An immediate consequence of the independentness is the following corollary that is first proved in
[39].
Corollary 4.3 If a node v i in V 2 is independent to all the nodes in V 1 , then s i.e., the values
of the ith row of C will not be modified in the partial LU factorization.
We now modify our threshold tolerance reordering strategy slightly to a diagonal threshold
strategy, similar to that discussed in [39]. We assume that the node v i is in V 1 if ja ii j - ffl and D
is still a diagonal matrix. With such a modification, we have jd m. Denote by
the size of the largest elements in absolute value of A.
Proposition 4.4 The size of the elements of the Schur complement matrix A 1 is bounded by M(1+
mM=ffl).
Proof. Starting from Equation (6):
+mM=ffl):Proposition 4.4 shows that the size of the elements of the Schur complement matrix cannot
grow uncontrollably if ffl is large enough. This result indicates that our first level (I)LU factorization
is stable.
As we hinted previously, we will be interested in recursively applying our strategy to the successive
Schur complement matrices. We may assume that the matrix A is presparsified so that small
nondiagonal elements are removed. To be more specific, for the parameter - used in the ILUT fac-
torization, we assume min 1-i;j-n fja ij jg - for all nonzero elements of A, except for possibly the
diagonal elements. With some additional assumptions, we can have a lower bound on the variation
of the diagonal values of the Schur complement matrix A 1 .
Proposition 4.5 Suppose ja ij j - for all nonzero elements of the matrix A, and suppose that either
c ik f ki =d k - 0 or c ik f ki =d k - 0 holds for all 1 - k - m. Then
card
where are the index sets of the nonzero elements of the ith row of the E submatrix
and the ith column of the F submatrix, respectively. card(V ) denotes the cardinality of a set V .
Proof. If either c ik f ki =d k - 0 or c ik f ki =d k - 0 holds for all 1 - k - m, we have
The kth term in the right-hand side sum of (7) is nonzero if and only if both e ik and f ki are nonzero.
This happens if and only if k 2 Nz(E
Note that jc ik j -; jf ki j - and jd k m. It follows that
card
:It is implicitly assumed that ffl ! M . In practice, ffl is small so that the set V 1 may be large
enough to avoid constructing a large Schur complement matrix. Denote
card
By the motivation of the diagonal threshold strategy, the value of jc ii j is zero or very small. Thus
the size of js ii j can be considered as being close to \Delta i .
Corollary 4.6 Under the conditions of the Proposition 4.5, if c
Corollary 4.6 shows that if the ith diagonal element of A 1 is zero in A and if the set Nz(E i
is nonempty, then the size of the ith diagonal element is nonzero in the Schur complement. Thus,
under these conditions, a zero pivot is removed. In fact, the cardinality of Nz(E seems to
be the key factor to remove zero diagonal elements.
It is difficult to derive more useful bounds for general sparse matrices. If certain conditions
are given to restrict the class of matrices under consideration, it is possible to obtain more realistic
bounds to characterize the size of the elements of the Schur complement matrix, especially the size
of its diagonal elements.
General submatrix D. For general submatrix D corresponding to the factorization (4), it is easy
to see that, if the jth column of the submatrix F is zero, the jth column of the submatrix L \Gamma1 F is
zero. Hence, Proposition 4.1 carries over to the general case.
At this moment, we are unable to show results analogous to Propositions 4.4 and 4.5 for general
submatrix D. However, it can be argued heuristically that, if D is not a diagonal matrix, the
cardinality of the set Nz(E i likely to be larger than that of Nz(E
Size of k(LU) \Gamma1 k. Let us consider the quality of preconditioning in a nonstandard way. Denote by
the error (residual) matrix of the ILU factorization. At each iteration, the preconditioning step solves
for -
w the system
where r is the residual of the current iterate. In a certain sense, we can consider -
w as an approximate
to the correction term of the current iterate. The quality of the preconditioning step (9) can be
judged by comparing (9) with the exact or perfect preconditioning step
If Equation (10) could be solved to yield the exact correction term w, the preconditioned iterative
method would converge in one step. Of course, solving the Equation (10) is as hard as solving the
original system (1). However, we can measure the relative difference in the correction term when
approximating the Equation (10) by the Equation (9). This difference may tell us how good the
preconditioning step (9) approximates the exact preconditioning step (10). The following proposition
is motivated by the work of Kershaw [23].
Proposition 4.7 Suppose the matrix A and the factor LU from the incomplete LU factorization are
nonsingular, the following inequality holds:
wk
kwk
for any consistent norm k \Delta k.
Proof. It is obvious that r 6= 0, otherwise the iteration would have converged. The nonsingularity
of A implies that w 6= 0. Note that -
Equation (10), we have
It follows that, for any consistent norm,
The desired result (11) follows immediately by dividing kwk on both sides. 2
It is well known that the size of the error matrix E directly affects the convergence rate of the
preconditioned iterative methods [16]. Proposition 4.7 shows that the quality of a preconditioning
step is directly related to the size of both (LU) \Gamma1 and R. A high quality preconditioner must be
accurate; i.e., it must have an error matrix that is small in size. A high quality preconditioner must
also have a stable factorization and stable triangular solves; i.e., the size of (LU) \Gamma1 must be small.
Since the condition estimate, condest is a lower bound for k(LU) \Gamma1 k1 , it should
provide some information about the quality of the preconditioner and may be used to measure the
stability of the LU factorization and of the triangular solves.
5 Multilevel Dual Reordering and ILU Factorization
Based on our previous analyses, the size of a diagonal element of the matrix A 1 is likely to be
larger than that of the same element in A. 2 We can apply Algorithm 2.1 to A 1 and repeat on
A 1 the procedures that were applied to A. This process may be repeated for a few times until
2 This is obviously false for an M-matrix. However, there will be no Schur complement matrix at all if A is an
all small diagonal elements are modified to large values, or until the last Schur complement matrix
is small enough that an ILU factorization with a complete pivoting strategy can be implemented
inexpensively. Since the number of small or zero pivots in the last Schur complement matrix is small,
a third strategy is to replace them by a large value. This will not introduce too much error to the
overall factorization. Given a maximum level L and denote A the multilevel dual reordering
strategy and ILU factorization can be formulated as Algorithm 5.1.
Algorithm 5.1 Multilevel dual reordering and ILU factorization.
1. Given the parameters -; p; ffl; L
2. For
3. Run Algorithm 2.1 with ffl to find permutation matrices
and P jg
4. Perform matrix permutation A
5. If no small pivot has been found, then
6. Apply ILUT(p; -) to A j and exit
7. Else
8. Apply a partial ILU factorization to A j
9. to yield a Schur complement matrix A j+1
10. End if
11. End do
12. Apply ILUTP or a stabilized ILUT to AL if AL exists
The ILU preconditioner constructed by Algorithm 5.1 is structurally similar to the BILUTM
preconditioner in [38]. The difference is that we do not construct a block independent set for the D j
submatrix. Instead, we set up a diagonal measure constraint and employ a graph reordering scheme
to reduce fill-in. The emphasis of this paper is on solving indefinite matrices by removing small pivots.
It can be seen, if L levels of reduction are performed, the resulting ILU preconditioner has the
following
The application of the preconditioner can be done by a level by level forward elimination, followed
by a level by level backward substitution. There are also permutations and inverse permutations to
be performed, specific procedures depend on implementations. For detailed descriptions, we refer to
[37, 38].
6 Numerical Experiments
Standard implementations of multilevel preconditioning methods have been described in detail in
[33, 37, 38]. We used full GMRES as the accelerator [35]. We tested three preconditioners: standard
ILUT of [32], a column pivoting variant ILUTP [32], and the multilevel dual reordering preconditioner
designed in this paper, abbreviated as MDRILU (multilevel dual reordering ILU factorization). All
preconditioners used a safeguard (stabilization) procedure by replacing a zero pivot with (0:0001+-)r i ,
where r i was computed as the average nonzero values of the row in question. They were used as right
preconditioners for GMRES [34]. The main parameters used in all three preconditioners are the pair
(p; -) in the double dropping strategy. ILUTP needs another parameter 0 - oe - 1 to control the
actual pivoting. A nondiagonal element a ij is a candidate for a permutation only when oeja
It is suggested that reasonable values of oe are between 0:5 and 0:01, with 0:5 being the best in
many cases [34, p. 295]. MDRILU also needs another parameter ffl to enforce the diagonal threshold
reordering as in Algorithm 5.1. The maximum possible level number in MDRILU was
levels of dual reorderings the Schur complement A 10 is not empty, a stabilized ILUT factorization
was employed to factor A 10 . 3
For all linear systems, the right-hand side was generated by assuming that the solution is a
vector of all ones. The initial guess was a vector of some random numbers. The iteration was
terminated when the 2-norm of the residual was reduced by a factor of 10 7 . We also set an upper
bound of 100 for the full GMRES iteration. A symbol"-" indicates lack of convergence.
In all tables with numerical results, "iter" shows the number of preconditioned GMRES iter-
ations; "spar" shows the sparsity ratio which is the ratio between the number of nonzero elements
of the preconditioner to that of the original matrix; "prec" shows the CPU time in seconds spent in
constructing the preconditioners; is the condition estimate of the preconditioners
as introduced in Section 1. Since these ILU preconditioners approach direct solvers as
robustness with respect to the memory cost (sparsity ratio). We remark
that our codes were not optimized and they computed and outputed information such as the number
of zero diagonals, smallest pivots, ets. Consequently, the CPU times reported in this paper only have
relative meaning. Note that the solution time at each iteration is mainly the cost of the matrix (both
A and the preconditioner) vector products and is thus proportional to the product of the iteration
count and the sparsity ratio, i.e., solution time
The numerical experiments were conducted on a Power-Challenge XL Silicon Graphics workstation
equipped with 512 MB of main memory, one 190 MHz R10000 processor, and 1 MB secondary
cache. We used Fortran 77 programming language in 64 bit arithmetic computation.
Test matrices. Three test matrices were selected from different applications. Table 1 contains
simple descriptions of the test matrices. They have been used in several other papers [6, 9, 39, 46].
None of the three matrices has a zero diagonal.
Matrix order nonzeros description
buckling problem for container model
simulation
WIGTO966 3 864 238 252 Euler equation model

Table

1: Simple descriptions of the test matrices.
3 We found stabilized ILUT was better than ILUTP for solving the last system. We did not implement an ILUT
factorization with a full pivoting strategy.
WIGTO966 matrices. The WIGTO966 matrix 4 was supplied by L. Wigton from Boeing Com-
pany. It is solvable by ILUT with large values of p [6]. This matrix was also used to compare BILUM
with ILUT in [36], and BILUTM with ILUT in [38], and to test point and block preconditioning
techniques in [8, 9]. Since ILUT requires very large amount of fill-in to converge, the WIGTO966
matrix is ideal to test alternative preconditioners and to show the least memory that is required for
convergence. For example, BILUM (with GMRES(10)) was shown to be 6 times faster than ILUT
with only one-third of the memory required by ILUT [36]. BILUTM (with GMRES(50)) converged
almost 5 times faster and used just about one-fifth of the memory required by ILUT [38]. Table 2 lists
results from several runs to compare MDRILU and ILUT. It shows that MDRILU could converge
with low sparsity ratios, as low as 0:94. The threshold parameter ffl was in a fixed range when the
other parameters p and - changed. For all the values of p and - tested in Table 2, ILUT did not
converge. We found that there was no very small pivot, the size of the smallest pivot in all tests in

Table

was 1.19e-5. But the condition estimates for ILUT were very large, the smallest condest value
is 1.1e+82, indicating unstable triangular solves had resulted during the factorization and solution
processes.
MDRILU ILUT
iter prec spar cond iter prec spar cond
50 1.0e-3 0.38 27 4.92 2.17 4.4e+4 - 9.9e+116
50 1.0e-4 0.38 25 7.48 2.55 2.7e+4 - 2.7e+91

Table

2: Comparison of MDRILU and ILUT for solving the WIGTO966 matrix.
We further compared ILUTP and ILUT and list the results in Table 3. We see that ILUTP
is more robust than ILUT for solving the WIGTO966 matrix. ILUT required high sparsity ratios
to converge. For those cases, ILUTP was able to converge with fewer iterations. When we chose
failed to converge, but ILUTP converged in 49 iteration with a sparsity
ratio 3:06. Notice that both ILUTP and ILUT did not converge with
MDRILU could converge with these parameters. We point out that the condition estimates of ILUTP
are much smaller than those of ILUT. This implies that ILUTP did stabilize the ILU factorization
process with a column pivoting strategy, although there was no very small pivot in the factorization.
The results of Table 3 also show that the additional cost of implementing ILUTP is not high in
this test. However, as far as solving the WIGTO966 matrix is concerned, computing an MDRILU
preconditioner is much cheaper than computing either an ILUT or an ILUTP preconditioner.
RAEFSKY4 matrices. The RAEFSKY4 matrix 5 was supplied by H. Simon from Lawrence
Berkeley National Laboratory (originally created by A. Raefsky from Centric Engineering). This is
4 The WIGTO966 matrix is available from the author.
5 The RAEFSKY4 matrix is available online from the University of Florida Sparse Matrix Collection at
http://www.cise.ufl.edu/~davis/sparse.
ILUTP ILUT
iter prec spar cond iter prec spar cond
100 1.0e-4 0.50 34 22.90 3.08 2.3e+5 - 3.0e+69
300 1.0e-3 0.10 9 44.98 7.39 1.2e+4 74 51.52 7.91 1.3e+8

Table

3: Comparison of ILUTP and ILUT for solving the WIGTO966 matrix.
probably the hardest one in the total of 6 RAEFSKY matrices. Figure 2 shows the convergence
history of three preconditioners with 1.0e-4. The other parameters were
for MDRILU and oe = 0:03 for ILUTP. We see that both ILUT and ILUTP did not have much
convergence in 100 iterations; MDRILU converged in 13 iterations.
2-norm
residual
iterations
RAEFSKY4 Matrix
DRILUDRILU (dashed line)
ILUTP (dashdot line)
ILUT (solid line)

Figure

2: Convergence history of preconditioned GMRES for solving the RAEFSKY4 matrix.
In

Figure

3 we plotted the iteration counts (left part) and the values of condition estimate
(right part) of the MDRILU preconditioner with different values of the threshold parameter ffl, keeping
We found that the iteration count and the condition estimate were linked
to each other. A large value of condition estimate is usually accompanied by a large iteration count
of MDRILU. We also see that the convergence rates of MDRILU are not very sensitive to the choice
of the value of ffl. For 0:38 - ffl - 0:78, MDRILU gave very similar performance.
UTM5940 matrix. The UTM5940 matrix 6 is the largest matrix from the TOKAMAK collection
and was provided by P. Brown of Lawrence Livermore National Laboratory. Table 4 contains a few
runs with MDRILU and ILUT with different sparsity ratios. It is clear that MDRILU is more efficient
than ILUT when the sparsity ratios are low. The results are also consistent with other test results,
6 The UTM5940 matrix is available from online the MatrixMarket of the National Institute of Standards and Tech-
RAEFSKY4 Matrix
iterations
epsilon
0.850015002500condest
epsilon
RAEFSKY4 Matrix

Figure

3: Iteration counts (left) and condition estimates (right) of MDRILU with different values of
ffl for solving the RAEFSKY4 matrix.
indicating that MDRILU is able to solve this problem with less storage cost than ILUT. If sufficient
memory space is available, ILUT may be efficient in certain cases. Note that if both MDRILU and
ILUT converge with similar iteration counts, MDRILU is more expensive to construct than ILUT.
MDRILU ILUT
iter prec spar cond iter prec spar cond
50 1.0e-4 0.30 42 7.49 6.26 2.2e+7 86 3.67 5.72 1.3e+7

Table

4: Comparison of MDRILU and ILUT for solving the UTM5940 matrix.

Figure

4 shows the convergence history of MDRILU with different values of dropping tolerance
- to solve the UTM5940 matrix, keeping We note that the number of
iterations did not change very much when - changed from 1.0e-2 to 1.0e-5 and the sparsity ratio
changed from 2:67 to 4:15. It seems that MDRILU worked quite well with a relatively strict dropping
tolerance.
FIDAP matrices. The FIDAP matrices 7 were extracted from the test problems provided in
the FIDAP package [21]. They were generated by I. Hasbani of Fluid Dynamics International and
B. Rackner of Minnesota Supercomputer Center. The matrices were resulted from modeling the incompressible
Navier-Stokes equations and were generated using whatever solution method was specified
in the input decks. However, if the penalty method was used, there is usually a corresponding
7 All FIDAP matrices are available online from the MatrixMarket of the National Institute of Standards and Tech-
2-norm
residual
iterations
Matrix
solid line:
dashed line:
dashdot line:
dotted line:

Figure

4: Convergence history of MDRILU with different values of dropping tolerance - for solving
the UTM5940 matrix.
FIDAPM matrix, which was constructed using a fully coupled solution method (mixed u-p formula-
tion). The penalty method gives very ill conditioned matrices, whereas the mixed u-p method gives
indefinite, larger systems (they include pressure variables).
Many of these matrices contain small or zero diagonal values. 8 The zero diagonals are due to
the incompressibility condition of the Navier-Stokes equations [9]. The substantial amount of zero
diagonals makes these matrices indefinite. It is remarked in [6] that the FIDAP matrices are difficult
to solve with ILU preconditioning techniques, which require high level of fill-in to be effective and the
performance of the preconditioners is unstable with respect to the amount of fill-in. Many of them
cannot be solved by the standard BILUM preconditioner and in some cases, even the construction of
BILUM failed due to the occurrence of very ill conditioned blocks. Nevertheless, some of them may
be solved by the enhanced version of BILUM using singular value decomposition based regularized
inverse technique and variable block size [40].
The details of all of the largest 31 FIDAP matrices (n ? 2000) are listed in Table 5 and the
corresponding test results are given in Table 6. The second column of Table 6 lists the number of
zero diagonals of the given matrix. In our tests, we first set
0:5; 0:3; 0:1; 0:01. If none of these ffl values showed any promise, we increased the p value or decreased
the - value. If for a given pair of (p; - ), MDRILU with a certain value of ffl converged or showed some
convergence, we adjusted the value of ffl to get improved convergence rates if possible. However, there
was no effort made to find the best parameters. We stopped refining the parameters when we found
the iteration count was reasonable and the sparsity ratio was not high, or the computations took too
much time in case of large matrices. Once MDRILU was tested, the same pair (p; -) was used to test
ILUTP and ILUT. For ILUTP, we varied the value of oe analogously to what we did to choose the
value of ffl.

Table

6 shows that MDRILU can solve 27 out of the 31 largest FIDAP matrices. To the best of
8 The FIDAP matrices have structural zeros added on the offdiagonals to make them structurally symmetric. Structural
zeros were also added to the diagonals.
Matrix order nonzeros description
developing flow in a vertical channel
impingment cooling
flow over multiple steps in a channel
flow in lid-driven wedge
FIDAP015 6 867 96 421 spin up of a liquid in an annulus
turbulent flow over a backward-facing step
developing pipe flow, turbulent
attenuation of a surface disturbance
coating
convection
two merging liquids with an interior interface
turbulent flow in axisymmetric U-bend
species deposition on a heated plate
FIDAP035 19 716 218 308 turbulent flow in a heated channel
FIDAP036 3 079 53 851 chemical vapor deposition
flow of plastic in a profile extrusion die
flow past a cylinder in free stream
natural convection in a square enclosure
developing flow in a vertical channel
impingment cooling
flow over multiple heat sources in a channel
FIDAPM11 22 294 623 554 3D steady flow, head exchanger
FIDAPM15 9 287 98 519 spin up of a liquid in an annulus
turbulent flow is axisymmetric U-bend
radiation heat transfer in a square cavity
FIDAPM37 9 152 765 944 flow of plastic in a profile extrusion die

Table

5: Description of the largest 31 FIDAP matrices.
our knowledge, this is the first time that so many FIDAP matrices were solved by a single iterative
technique. (20 were solved in [40], in [46], 9 in [39], and 8 in [9].) In Table 6 the term "unstable"
means that convergence was not reached in 100 iterations and the condition estimate was greater than
Similarly the term "inaccurate" means that convergence was not reached, but the condition
estimate did not exceed 10 15 . They are categorized according to Chow and Saad's arguments [9]. We
remark that the results of "inaccurate" or "unstable" in Table 6 do not indicate that ILUT or ILUTP
can or cannot solve the given matrices with different parameters. The results only mean that they
did not converge with the parameters that made MDRILU converge. It is worth pointing out that,
in several tests, we observed that ILUTP encountered zero pivots when ILUT did not.
Although we allowed 10 levels of maximum dual reorderings to be performed, there were very
few cases that 10 levels of reorderings were actually needed. In most cases, 3 to 4 levels of dual
reorderings were performed for the FIDAP matrices. In many cases, the first Schur complement
MDRILU ILUTP ILUT
Matrix zero-d p - ffl iter spar iter spar iter spar
unstable unstable
unstable unstable
FIDAP026 457 20 1.0e-4 0.30 84 0.77 unstable unstable
FIDAP036 504 20 1.0e-4 0.10 23 1.75 83 1.91 unstable
FIDAPM07 432 300 1.0e-4 0.20 78 6.66 80 7.71 inaccurate
FIDAPM08 780 20 1.0e-4 0.10 25 1.70 78 2.22 unstable
unstable unstable
43 7.51 21 7.37 unstable
28 14.38 11 7.47 13 7.46
43 3.61 unstable unstable

Table

Solving the FIDAP matrices by MDRILU, ILUTP and ILUT.
matrix did not have any zero diagonal, even if the original matrix A did have many zero diagonals.
We listed in Table 7 those matrices that did have zero diagonals in their Schur complement matrices.
For all the FIDAP matrices solved by MDRILU, only the FIDAP026 matrix had 12 zero diagonals
in the last Schur complement A 5 . The test results show that the multilevel dual reordering strategy
does have the effect of removing small and zero pivots from ILU factorizations.
Remarks. Ironically, the four matrices, FIDAP011, FIDAP015, FIDAP018, and FIDAP035, that
were not solved by MDRILU do not have any zero diagonals. They may be solved by ILUT with
small values of - . Some of them may even be solved by GMRES without preconditioning if enough
iterations are allowed. We think this is because these matrices are very nonsymmetric and the
preconditioned matrices were worse conditioned than the original matrices, causing GMRES iteration
to converge extremely slowly. One of our strong feeling in these numerical experiments is that, in
general, MDRILU does not seem to work well when - is very small. Large values of p usually improve
convergence. This observation can be seen in Figure 5 which depicts the convergence history of
Matrix A 0 A 1 A 3 A 4 A 5

Table

7: Number of zero diagonals in the Schur complement matrices.
MDRILU for solving the largest FIDAP matrix, FIDAPM11. We used tested
two values of 1.0e-3. It is clear that more accurate (in terms of dropping tolerance)
ILU factorization does not help and sometimes hampers convergence. Good values for the parameter
ffl are between 0:1 and 0:5. For most problems, the performance of MDRILU is not very sensitive to
the choice of ffl, as long as it is in the range of 0:1 and 0:5.
2-norm
residual
iterations
FIDAPM11 Matrix
solid line:
dashed line:

Figure

5: Convergence history of MDRILU with different values of dropping tolerance - for solving
the FIDAPM11 matrix.
7 Conclusion
We have proposed a multilevel dual reordering strategy for constructing robust ILU preconditioners
for solving general sparse indefinite matrices. This reordering strategy is combined with a partial
ILU factorization procedure to construct recursive Schur complement matrices. The preconditioner
is a multilevel ILU preconditioner. However, the constructed preconditioner (MDRILU) is different
from all existing multilevel preconditioners in a fundamental concept [37, 47]. MDRILU never intends
to utilize any traditional multilevel property, it uses the Schur complement approach solely for the
purpose of removing small pivots.
We conducted analyses on simplified model problems to find out how the size of the small diagonal
elements and other elements is modified when these elements become the elements of the Schur
complement matrix. We gave an upper bound for the size of general elements of the Schur complement
matrix to show that their size will not grow uncontrollably if a suitable threshold reordering based
on the diagonal dominance measure is implemented. We also showed that under certain conditions,
a zero or very small diagonal element is likely to be modified to favor a stable ILU factorization by
the Schur complement procedure.
We further studied the quality of a preconditioning step. We showed that the quality of a
preconditioning step is directly related to the size of both (LU) \Gamma1 and R (the error matrix). Hence,
a high quality preconditioner must have a stable ILU factorization and stable triangular solves, as
well as a small size error matrix. In other words, both accuracy and stability affect the quality of a
preconditioner.
We performed numerical experiments to compare MDRILU with two popular ILU precondi-
tioners. Our numerical results show that MDRILU is much more robust than both ILUT and ILUTP
for solving most indefinite matrices under current consideration. The most valuable advantage of
MDRILU is that it can construct a sparse high quality preconditioner with low storage cost. The
preconditioners computed by MDRILU are more stable than those computed by ILUT and ILUTP,
thanks to the ability of MDRILU to remove (not replace) the small diagonal values.
Both analytic and numerical results strongly support our conclusion that the multilevel dual
reordering strategy developed in this paper is a very useful strategy to construct robust ILU preconditioners
for solving general sparse indefinite matrices. Due to the time and space limit, we have not
tested other graph reordering algorithms in the multilevel dual reordering algorithm. Some of the
popular reordering strategies such as Cuthill-McKee and reverse Cuthill-McKee algorithms may be
useful in such applications to further improve quality of the ILU preconditioner. However, we fell the
robustness of MDRILU is mainly a result of using threshold tolerance reordering strategy and partial
ILU factorization to remove small pivots. The difference arising from using different graph algorithm
may be significant in terms of the number of iterations. But such a difference is unlikely to alter the
stability problem in a systematic manner in the ILU factorization.



--R

Comparison of fast iterative methods for symmetric systems.
Incomplete factorization methods for fully implicit simulation of enhanced oil recovery.
Orderings for incomplete factorization preconditioning of nonsymmetric problems.
An incomplete-factorization preconditioning using red-black ordering
Parallel elliptic preconditioners: Fourier analysis and performance on the Connection machine.

On preconditioned Krylov subspace methods for discrete convection-diffusion problems
An object-oriented framework for block preconditioning
Experimental study of ILU preconditioners for indefinite matrices.
Weighted graph based ordering techniques for preconditioned conjugate gradient methods.
Ordering methods for preconditioned conjugate gradient methods applied to unstructured grid problems.
SOR as a preconditioner.
Stability and spectral properties of incomplete factorization.
On parallelism and convergence of incomplete LU factorizations.
Direct Methods for Sparse Matrices.
The effect of reordering on preconditioned conjugate gradients.
The effect of reordering on the preconditioned GMRES algorithm for solving the compressible Navier-Stokes equations
A stability analysis of incomplete LU factorization.
Relaxed and stabilized incomplete factorization for nonselfadjoint linear systems.
Ordering techniques for the preconditioned conjugate gradient method on parallel computers.
FIDAP: Examples Manual
Computer Solution of Large Sparse Positive Definite Systems.
On the problem of unstable pivots in the incomplete LU-conjugate gradient method
Conjugate gradient methods and ILU preconditioning of non-symmetric matrix systems with arbitrary sparsity patterns
Comparative analysis of the Cuthill-McKee and the reverse Cuthill-McKee ordering algorithms for sparse matrices
Ordering strategies for modified block incomplete factorizations.
An incomplete factorization technique for positive definite linear systems.
An iterative solution method for linear systems of which the coefficient matrix is a symmetric M-matrix
On the stability of the incomplete LU-factorization and characterizations of H-matrices
Ordering methods for approximate factorization preconditioning.
Orderings for conjugate gradient preconditionings.
ILUT: a dual threshold incomplete LU preconditioner.
ILUM: a multi-elimination ILU preconditioner for general sparse matrices
Iterative Methods for Sparse Linear Systems.
GMRES: a generalized minimal residual method for solving non-symmetric linear systems
Domain decomposition and multi-level type techniques for general sparse linear systems
BILUM: block versions of multielimination and multilevel ILU preconditioner for general sparse linear systems.
BILUTM: a domain-based multi-level block ILUT preconditioner for general sparse matrices
Diagonal threshold techniques in robust multi-level ILU preconditioners for general sparse linear systems
Enhanced multi-level block ILU preconditioning strategies for general sparse linear systems
A multi-level preconditioner with applications to the numerical simulation of coating problems
On the stability of the incomplete Cholesky decomposition for a singular perturbed problem
Incomplete LU preconditioners for conjugate-gradient-type iterative methods
Iterative solution methods for certain sparse linear systems with a non-symmetric matrix arising from PDE-problems
Stabilized incomplete LU-decompositions as preconditionings for the Tchebycheff iteration
Preconditioned Krylov subspace methods for solving nonsymmetric matrices from CFD applications.
A grid based multilevel incomplete LU factorization preconditioning technique for general sparse matrices.
--TR

--CTR
Wang , Jun Zhang, A new stabilization strategy for incomplete LU preconditioning of indefinite matrices, Applied Mathematics and Computation, v.144 n.1, p.75-87, 20 November
Kai Wang , Jun Zhang , Chi Shen, Parallel Multilevel Sparse Approximate Inverse Preconditioners in Large Sparse Matrix Computations, Proceedings of the ACM/IEEE conference on Supercomputing, p.1, November 15-21,
Jeonghwa Lee , Jun Zhang , Cai-Cheng Lu, Incomplete LU preconditioning for large scale dense complex linear systems from electromagnetic wave scattering problems, Journal of Computational Physics, v.185 n.1, p.158-175, 10 February
Chi Shen , Jun Zhang, Parallel two level block ILU Preconditioning techniques for solving large sparse linear systems, Parallel Computing, v.28 n.10, p.1451-1475, October 2002
Chi Shen , Jun Zhang , Kai Wang, Distributed block independent set algorithms and parallel multilevel ILU preconditioners, Journal of Parallel and Distributed Computing, v.65 n.3, p.331-346, March 2005
Michele Benzi, Preconditioning techniques for large linear systems: a survey, Journal of Computational Physics, v.182 n.2, p.418-477, November 2002
sparse matrices;multilevel incomplete LU preconditioner;incomplete LU factorization;reordering strategies