Ordering, Anisotropy, and Factored Sparse Approximate Inverses.
We consider ordering techniques to improve the performance of factored sparse approximate inverse preconditioners, concentrating on the AINV technique of M. Benzi and M. T\r{u}ma. Several practical existing unweighted orderings are considered along with a new algorithm, minimum inverse penalty (MIP), that we propose. We show how good orderings such as these can improve the speed of preconditioner computation dramatically and also demonstrate a fast and fairly reliable way of testing how good an ordering is in this respect. Our test results also show that these orderings generally improve convergence of Krylov subspace solvers but may have difficulties particularly for anisotropic problems. We then argue that weighted orderings, which take into account the numerical values in the matrix, will be necessary for such systems. After developing a simple heuristic for dealing with anisotropy we propose several practical algorithms to implement it. While these show promise, we conclude a better heuristic is required for robustness.
Introduction
. Consider solving the system of linear equations:
where A is a sparse nn matrix. Depending on the size of A and the nature of the
computing environment an iterative method, with some form of preconditioning to
speed convergence, is a popular choice. Approximate inverse preconditioners, whose
application requires only (easily parallelized) matrix-vector multiplication, are of particular
interest today. Several methods of constructing approximate inverses have been
proposed (e.g., [2, 3, 9, 20, 22, 24]), falling into two categories: those that directly
form an approximation to A -1 and those that form approximations to the inverses of
the matrix's LU factors. This second category currently shows more promise than the
first for three reasons. First, it is easy to ensure that the factored preconditioner is
nonsingular, simply by making sure both factors have nonzero diagonals. Second, the
factorization appears to allow more information per nonzero to be stored, improving
convergence [4, 8]. Third, the setup costs for creating preconditioners can often be
much less [4].
However, unlike A -1 itself, the inverse LU factors are critically dependent on the
ordering of the rows and columns-indeed, they will not exist in general for some
orderings. Even in the case of an SPD matrix, direct methods have shown how important
ordering can be. Thus any factored approximate inverse scheme must handle
# Received by the editors March 6, 1998; accepted for publication (in revised form) March 18, 1999;
published electronically November 17, 1999. This research was supported by the Natural Sciences and
Engineering Council of Canada, and Communications and Information Technology Ontario (CITO),
funded by the Province of Ontario.
http://www.siam.org/journals/sisc/21-3/33584.html
Stanford University, SCCM program, Room 284, Gates Building 2B, Stanford, CA 94305
# Department of Computer Science, University of Waterloo, Waterloo, ON, N2L 3G1, Canada
(wptang@riacs.edu).
868 ROBERT BRIDSON AND WEI-PAI TANG
ordering with thought. In particular, for an e#ective preconditioner an ordering that
minimizes the size of the "dropped" entries is needed-decreasing the error between
the approximate inverse factors and the true ones (see [14] for a discussion of this in
the context of ILU).
In this paper we focus our attention on the AINV algorithm [3], which, via implicit
Gaussian elimination with small-element dropping, constructs a factored approximate
where Z and W are unit upper triangular and D -1 is diagonal. However, the purely
structural results presented in section 2 apply equally to other factored approximate
inverse schemes. Whether the numerical results carry over is still to be determined.
For example, conflicting evidence has been presented in [5] and [16] about the e#ect
on FSAI [22], which perhaps will be resolved only when the issue of sparsity pattern
selection for FSAI has been settled.
Some preliminary work in studying the e#ect of ordering on the performance of
AINV has shown promising results [3]. more recent work by the same authors is
[5].) We carry this research forward in sections 2 and 3, realizing significant improvements
in the speed of preconditioner computation and observing some beneficial e#ects
on convergence but noting that structural information alone is not always enough. We
then turn our attention to anisotropic problems. For the ILU class of preconditioners
it has been determined that orderings which take the numerical values of the matrix
into account are useful-even necessary (e.g., [10, 11, 12, 14]). Sections 4-6 try to answer
the question of whether the same thing holds for factored approximate inverses.
The appendix contains the details of our test results.
2. Unweighted orderings. Intuitively, the smaller the size of the dropped portion
from the true inverse factors, the better the approximate inverse will be. We will
for now assume that the magnitudes of the inverse factors' nonzeros are distributed
roughly the same way under di#erent orderings. (Our experience shows this is a fairly
good assumption for typical isotropic problems, but as we shall see later, this breaks
down for anisotropic matrices in particular.) Then we can consider the simpler problem
of reducing the number of dropped nonzeros, instead of their size. Of course, for
sparsity we also want to retain as few nonzeros as possible; thus we really want to
reduce the number of nonzeros in the exact inverse factors-a quantity we call inverse
factor fill, or I F fill.
Definition 2.1. Let A be a square matrix with a triangular factorization
LU . The I F fill of A is defined to be the total number of nonzeros in the inverses of
L and U, assuming no cancellation in the forming of those inverses.
For simplicity we restrict our discussion to the SPD case, first examining I F
fill and then considering several existing ordering algorithms that may be helpful.
We finish the section by proposing a new ordering scheme, which we call MIP. The
application to the unsymmetric case is straightforward.
The following discussion makes use of some concepts from graph theory. The
graph of an n  n matrix A is a directed graph on n nodes labeled 1, . , n, with
an arc i # j if and only if A ij #= 0. A directed path, or dipath, is an ordered set
such that the arcs all exist-often this is
written as u 1 # u k . See chapter 3 of [17], for example, for further explanation.
From Gilbert [19] and Liu [23], we have the following graph theoretic characterization
of the structure of the inverse Cholesky factor.
ORDERING, ANISOTROPY, AND APPROXIMATE INVERSES 869
Theorem 2.2. Let A be an SPD matrix with Cholesky factor L. Then assuming
no cancellation, the structure of -T corresponds to the transitive closure 1 of
the graph of L T , that is, for i < j, we have Z ij #= 0 if and only if there is a dipath
from i to j in the graph of L T . Furthermore, this is the same structure as given by
the transitive closure of the elimination tree of A (the transitive reduction 2 of L T ).
Notice that the last structure characterization simply means that for i > j,
is an ancestor of i in the elimination tree. This allows
us to significantly speed the computation of the preconditioner given a bushy elimination
tree, as well as allowing for parallelism-for the calculation of column j in the
factors, only the ancestor columns need be considered (a coarser-grain version of this
parallelism via graph partitioning has been successfully implemented in [6] for exam-
ple). Another product of this characterization is a simple way of computing the I F
fill of an SPD matrix, obtained by summing the number of nonzeros in each column
of the inverse factor and multiplying by two for the other (transposed) factor.
Theorem 2.3. The I F fill of an SPD matrix is simply twice the sum of the
depths of all nodes in its elimination tree. In particular, the number of nonzeros in
column j of the true inverse factor Z is given by the number of nodes in the subtree
of the elimination tree rooted at j.
These results suggest orderings that avoid long dipaths in L T (i.e., paths in L T
with monotonically increasing node indices), as these cause lots of I F fill, quadratic
in their length. Alternatively, we are trying to get short and bushy elimination trees.
Another useful characterization of I F fill using notions from [17, 18] allows us
to do a cheap "inverse symbolic factorization"-determining the nonzero structure
of the inverse factors-without using the elimination tree, which is essential for our
minimum inverse penalty (MIP) ordering algorithm presented later.
Theorem 2.4. Z ij #= 0 if and only if i is reachable from j strictly through nodes
eliminated previous to i-or in terms of the quotient graph model, if i is contained in
a supernode adjacent to j at the moment when j is eliminated.
Based on the heuristic and results above, we now examine several existing orderings
which might do well and propose a new scheme to directly implement the
heuristic of reducing I F fill.
Red-black. The simplest ordering we consider is (generalized) red-black, where a
maximal independent set of ("red") nodes is ordered first and the remaining ("black")
nodes are ordered next according to their original sequence. In that initial red block,
there are no nontrivial dipaths and hence no o#-diagonal entries in Z.
Minimum degree. Benzi and Tuma have observed that minimum degree is
generally beneficial for AINV [3]. This is justified by noting that minimum degree
typically substantially reduces the height of the elimination tree and hence should
reduce I F fill.
As an aside, notice that direct-method fill-reducing orderings do not necessarily
reduce I F fill. For example, a good envelope ordering will likely give rise to a very
tall, narrow elimination tree-typically just a path-and thus give full inverse factors.
Again, it should be noted that this isn't necessarily a bad thing if the inverse factors
still have very small entries, but without using numerical information from the matrix
our experience is that envelope orderings do not manage this. This may seem at
1 The transitive closure of a directed graph G is a graph G # on the same vertices with an arc
vertices u and v that were connected by a dipath u # v in G.
2 The transitive reduction of a directed graph G is a graph G # with the minimum number of arcs
but still possessing the same transitive closure as G.
variance with the result in [13] (elaborated in [5] for factored inverses) that for banded
SPD matrices the rate of decay of the entries in the inverse has an upper bound that
decreases as the bandwidth decreases. However, decay was measured there in terms of
distance from the diagonal, which is really suitable only for small bandwidth orderings;
the results presented in section 5 of [13], measuring decay in terms of the unweighted
graph distance, should make theoretical progress possible. We will return to this issue
in sections 4-6.
Nested dissection (ND). On the other hand, by ordering vertex separators last
ND avoids any long monotonic dipaths and hence a lot of I F fill. Alternatively, in
trying to balance the elimination tree it reduces the sum of depths.
Minimum inverse penalty (MIP). Above we noted that we can very cheaply
compute the number of nonzeros in each column of Z within a symbolic factoriza-
tion. This allows to propose a new ordering, MIP, an analogue to minimum degree.
Minimum degree is built around a symbolic Cholesky factorization of the matrix, at
each step selecting the node(s) of minimum penalty to eliminate. The penalty was
originally taken to be the degree of the node in the partially eliminated graph; later
algorithms have used other related quantities including the external degree and approximate
upper bounds. In MIP we follow the same greedy strategy but we compute
a penalty for node v based on Zdeg v , the number of nonzeros in the column of Z
were v to be ordered next-the degree in the inverse Cholesky factor instead of the
well as on Udeg v , the number of uneliminated neighbors of node
v at the current stage of factorization (not counting supernodes). In our experiments
we found the function Penalty to be fairly e#ective. Further
research into a better penalty function is needed. Also, ideas from minimum degree
such as multiple elimination, element absorption, etc. might be suitable here.
3. Testing unweighted orderings. We used the symmetric part of the matrix
for all the orderings. Red-black was implemented with the straightforward greedy
algorithm to select a maximal independent set. Our minimum degree algorithm was
AMDBAR, a top-notch variant due to Amestoy, Davis, and Du# [1]. We wrote our
own ND algorithm that constructs vertex separators from edge separators given by a
multilevel bisection algorithm. This algorithm coarsens the graph with degree-1 node
compression and heavy-edge matchings until there are less than 100 nodes, bisects
the small graph spectrally according to the Fiedler vector [15], and uses a greedy
boundary-layer sweep to smooth in projecting back to the original. See [21], for
example, for more details on this point.
The appendix provides further details about our testing. The tables contain data
for both the unweighted orderings above and their weighted counterparts presented
below-ignore the lower numbers for now. In brief, we selected several matrices from
the Harwell-Boeing collection and tested them all with each ordering scheme. Table
4 gives the true I F fill for each matrix (or its symmetric part) and ordering. Tables
5-7 give the preconditioner performance. As the number of nonzeros allowed in the
preconditioner can have a significant e#ect on results, we standardized all our test
runs: in each box the left number is a report for when the preconditioner had as
many nonzeros as the matrix and the right number for when the preconditioner had
twice as many nonzeros.
In terms of I F fill reduction, AMDBAR, ND, and MIP are always the best three
by a considerable factor. ND wins 15 times, AMDBAR 5 times, and MIP 3 times,
with one tie. Red-black beats the natural ordering but not dramatically.
It is clear that ordering can help immensely for accelerating the computation of
ORDERING, ANISOTROPY, AND APPROXIMATE INVERSES 871
Inverse fill
Time
to
compute
preconditioner
Fig. 1. Correlation of I F fill and preconditioner computing time, normalized with respect to
the given (original) ordering.

Table
Average decrease in number of iterations over the test set, in percentages of the iterations taken
by the given ordering.
Fill Red-black AMDBAR ND MIP
the preconditioner. ND is the winner, followed by AMDBAR, MIP, and then red-black
quite a bit behind. The preconditioner computing time is closely correlated to
I F fill-see Figure 1. Thus calculating I F fill provides a fast and reasonably good test
to indicate how e#cient an (unweighted) ordering is for preconditioner computation-
perhaps not an important point if the iteration time dominates the setup time, but
this may be useful for applications where the reverse is true.
The e#ect of ordering on speed of solution is less obvious. The poor behavior of
PORES2, 3 SHERMAN2, and WATSON5 indicate that AINV probably isn't appropriate
(although if we had properly treated SHERMAN2 as a block matrix instead
it might have gone better). Notice in particular that sometimes lowering the drop
tolerance, increasing the size of the preconditioner and hopefully making it more
accurate, actually degrades convergence for these indefinite problems. From the remaining
matrices, we compared the average decrease in number of iterations over the
given ordering-see Table 1. Particularly given its problems with SAYLR4, WAT-
SON4, WATT1, and WATT2 red-black cannot be viewed as a good ordering. MIP is
overall the best, although it had a problem with WATT1, whereas the close contender
AMDBAR did fairly well on all but the di#cult three mentioned above. The good I F
fill reducing orderings all made worthwhile improvements in convergence rates, but it
is surprising that ND did the least considering its superior I F fill reduction. Clearly
our intuition that having fewer nonzeros to drop makes a better preconditioner has
merit, but it does not tell the entire story.
3 In [3] somewhat better convergence was achieved for PORES2, presumably due to implementation
di#erences in the preconditioner or its application.
872 ROBERT BRIDSON AND WEI-PAI TANG
4. Anisotropy. In the preceding test results we find several exceptions to the
general rule that AMDBAR, ND, and MIP perform similarly, even ignoring PORES2,
SHERMAN2, and WATSON5. For example, there are considerable variations for each
of ALE3D, BCSSTK14, NASA1824, ORSREG1, SAYLR4, and WATSON4. More
importantly, these variations are not correlated with I F fill; some other factor is at
work. Noticing that each of these matrices are quite anisotropic and recalling the
problems anisotropy poses for ILU, we are led to investigate weighted orderings.
We first develop a heuristic for handling anisotropic matrices. The goal in mind
is to order using the di#ering strengths of connections to reduce the magnitude of
the inverse factor entries. Then even if we end up with more I F fill (and hence drop
more nonzeros), the magnitude of the discarded portion of the inverse factors may be
smaller and give a more accurate preconditioner.
Again, we only look at SPD A. Let A = LDL T be the modified Cholesky fac-
torization, where L is unit lower triangular and D is diagonal. Then -T , and
since is zero on and below the diagonal hence nilpotent, we have
Then for
The nonzero entries in this sum correspond to the monotonically increasing di-
Our orderings should therefore avoid having
many such dipaths which involve large entries of L, as each one could substantially
increase the magnitude of Z's entries. Thus we want to move the large entries away
from the diagonal, so they cannot appear in many monotonic dipaths. In other words,
after a node has been ordered, we want to order so that its remaining neighbors with
strong L-connections come as late as possible after it.
For the purposes of our ordering heuristic, we want an easy approximation to L
independent of the eventual ordering chosen-something that can capture the order of
magnitude of entries in L but doesn't require us to decide the ordering ahead of time.
Assuming A has an adequately dominant diagonal without too much variation, we
can take the absolute value of the lower triangular part of A, symmetrically rescaled
to have a unit diagonal. (This can be thought of as a scaled Gauss-Seidel approxima-
tion.) Our general heuristic is then to delay strong connections in this approximating
matrix M defined by
An alternative justification of this heuristic, simply in the context of reducing the
magnitude of entries in L, is presented in [10].
Now consider a simple demonstration problem to determine whether this heuristic
could help. The matrix SINGLEANISO comes from a 5-point finite-di#erence
discretization on a regular 31  31 grid of the following PDE:
Here the edges of A (and M) corresponding to the y-direction are 1000 times heavier
than those corresponding to the x-direction. We try comparing two I F fill reducing
orderings. The first ordering ("strong-first") block-orders the grid columns with
ORDERING, ANISOTROPY, AND APPROXIMATE INVERSES 873
Strong-first ordering Weak-first ordering
Fig. 2. The two orderings of SINGLEANISO, depicted on the domain. Lighter shaded boxes
indicate nodes ordered later.

Table
Performance of strong-first ordering versus weak-first ordering for SINGLEANISO.
Ordering
Time to compute
preconditioner
Number of
iterations
Time for
iterations
Strong-first
Weak-first 0.38 25 0.25
nested dissection and then internally orders each block with nested dissection-this
brings the strong connections close to the diagonal. The second ("weak-first") block-
orders the grid rows instead, pushing the strong connections away from the diagonal,
delaying them until the last. These are illustrated in Figure 2, where each square of
the grid is shaded according to its place in the ordering.
Both orderings produce a reasonable I F fill of 103,682, with isomorphic elimination
trees. However, they give very di#erent performance at the first level of fill-see

Table

2. In all respects the weak-first ordering is significantly better than the strong-
first one.
In

Figure

3 we plot the decay of the entries in the inverse factors resulting from
the two orderings and show parts of those factors. The much smaller entries from the
weak-first ordering confirm our heuristic.
5. Weighted orderings. In our experience it appears that I F fill reduction
typically helps to also reduce the magnitude of entries in the inverse factors but blind
as it is to the numerical values in the matrix it can make mistakes such as allowing
strong connections close to the diagonal. In creating algorithms for ordering general
matrices, we thus have tried to simply modify the unweighted algorithms to consider
the numerical values.
Weighted nested dissection (WND). Consider the spectral bipartition al-
gorithm. Finding the Fiedler vector, the eigenvector of the Laplacian of the graph
with second smallest eigenvalue (see [15]), is equivalent to minimizing (over a space
orthogonal to the constant vectors):
Entries sorted by magnitude
Magnitudes
Comparing decay in factors
Strong-first
Weak-first
Close-up of strong-first
Close-up of weak-first
Fig. 3. Comparison of the magnitude of entries in inverse factors for the di#erent orderings of
SINGLEANISO. The close-up images of the actual factors are shaded according to the magnitude
of the nonzeros-darker means bigger.
(i,j) is an edge
We then make the bipartition depending on which side of the median each entry
lies. Notice that the closer together two entries are in value-i.e., the smaller
is-the more likely those nodes will be ordered on the same side of the cut. We would
like weakly connected nodes (where M ij is small) to be in the same part and the
strong connections to be in the edge cut, so we try minimizing the following weighted
quadratic sum:
(i,j) is an edge
where M is the scaled matrix mentioned above. This corresponds to finding the
eigenvector with second smallest eigenvalue of the weighted Laplacian matrix for the
graph defined by
(i, is an edge if and only if M ij #= 0, weight(i,
Thus we modify ND simply by changing the Laplacian used in the bipartition
step to this weighted Laplacian. Fortunately, our multilevel approach with heavy-
edge matchings typically will eliminate the largest o#-diagonal entries, as well as
substantially decreasing the size of the eigenproblem, making it easy to solve, so our
WND is very reasonable to compute.
OutIn. Our other weighted orderings are based on an intuition from finite-
di#erence matrices. We expect that the nodes most involved in long, heavy paths are
those near the weighted center of the graph (those with minimum weighted eccen-
tricity). The nodes least involved in paths are intuitively the ones on the weighted
ORDERING, ANISOTROPY, AND APPROXIMATE INVERSES 875
periphery of the graph. Thus OutIn orders the periphery first and proceeds to an
approximate weighted center ("from the outside to the inside"). To e#ciently find an
approximate weighted center we use an iterative algorithm.
Algorithm 1 (approximate weighted center).
. Do
- Calculate the M-weighted shortest paths and M-weighted distances to
other nodes from v i (where the distance is the minimum sum of weights,
given by M , along a connecting path).
- Select a node e i of maximum distance from v i .
Travel r of the way along a shortest path from v i to e i , saving the resulting
node as v i+1 .
. End loop when v as the approximate center.
Then our OutIn ordering is the following.
Algorithm 2 (OutIn).
. Compute
. Get an approximate weighted center c for M .
. Calculate the distances and shortest paths from all other nodes in M to c.
. Return the nodes in sorted order, with most distant first and c last.
Despite our earlier remark that envelope orderings might not be useful, the weight
information actually lets OutIn perform significantly better than the natural ordering
by reducing the size of the nonzeros in the factors if not the number of them. However,
why not combine OutIn with an I F fill reducing ordering to try to get the best of
both worlds? We thus test OutIn as a preprocessing stage before applying red-black,
minimum degree, or MIP. We note that the use of hash tables and other methods to
accelerate the latter two means that it's not true that the precedence set by OutIn
will always be followed in breaking ties for minimum penalty.
As an aside, we also considered modifying minimum degree and MIP with tie-breaking
directly based on the weight of a candidate node's connections to previously
selected nodes. Weighted tie-breaking (with RCM) has proved useful before in the
context of ILU [11]. However, for the significant extra cost incurred by this tie-
breaking, this achieved little here-it appears that a more global view of weights is
required when doing approximate inverses.
Before proceeding to our large test set, we verify that these orderings are behaving
as expected with another demonstration matrix. ANISO is a similar problem to
SINGLEANISO but with four abutting regions of anisotropy with di#ering directions
[12]-see

Figure

4 for a diagram showing the directions in the domain. As shown in

Table

3, the results for WND over ND and OutIn/MIP over MIP didn't change, but
there was a significant improvement in the other orderings.
6. Testing weighted results. We repeated the tests for our weighted orderings,
with results given in Tables 4-7. For unsymmetric matrices, we used |A|
define M in WND and OutIn, avoiding the issue of directed edges as with unweighted
orderings. In each box of the tables, the lower numbers correspond to the weighted
orderings; we have grouped them with the corresponding unweighted orderings for
comparison.
876 ROBERT BRIDSON AND WEI-PAI TANG
Fig. 4. Schematic showing the domain of ANISO. The arrows indicate the direction of the
strong connections.

Table
Performance of weighted orderings versus unweighted orderings for ANISO.
Ordering I F fill
Time to compute
preconditioner
Number of
iterations
Time for
iterations
Given 462 0.41 118 1.2
OutIn 266 0.31 77 0.76
Red-black 239 0.28 107 1.13
OutIn/RB 208 0.28 61 0.61
OutIn/AMD 69 0.15 48 0.5
WND 84 0.15
Only ND su#ered in preconditioner computing time-our spectral weighting appears
to be too severe, creating too much I F fill. However, it is important to note
that the increase in time is much less than that suggested by I F fill-indeed, although
WND gave several times more I F fill for ADD32, MEMPLUS, SAYLR4, SHERMAN4,
and WANG1, it actually allowed for slightly faster preconditioner computation. This
verifies the merit of our heuristic. Both the natural ordering and red-black benefited
substantially from OutIn in terms of preconditioner computation, and AMDBAR and
MIP didn't seem to be a#ected very much-this could quite well be a result of the
data structure algorithms which do not necessarily preserve the initial precedence set
by OutIn.
In terms of improving convergence, we didn't fix the problems with PORES2,
SHERMAN2, and WATSON5. These matrices have very weak diagonals anyhow, so
our heuristics probably don't apply. OutIn and OutIn/RB are a definite improvement
on the natural ordering and red-black, apart from on BCSSTK14 and WATSON4. The
e#ect of OutIn on AMDBAR and MIP is not clear; usually there's little e#ect, and
on some matrices (e.g., ALE3D and SAYLR4) it has an opposite e#ect on the two.
WND shows more promise, improving convergence over ND considerably for ALE3D,
BCSSTK14, ORSIRR1, ORSREG1, and SAYLR4. Its much poorer I F fill reduction
(generally by a factor of 4) gave it problems on a few matrices though.
ORDERING, ANISOTROPY, AND APPROXIMATE INVERSES 877
7. Conclusions. It is clear that I F fill reduction is crucial to the speed of preconditioner
computation, often making an order of magnitude di#erence. We also saw
that the I F fill of a matrix can be computed very cheaply and gives a good indication
of the preconditioner computation time, for unweighted orderings at least.
Reducing I F fill typically also gives a more e#ective preconditioner, accelerating
convergence-not only are the number of nonzeros in the true inverse factors de-
creased, but the magnitude of the portion that is dropped by AINV is reduced too.
However, although ND gave the best I F fill reduction, MIP gave the best acceleration
so care must be taken. It would be interesting to determine why this is so. Probably
several steps of ND followed by MIP or a minimum degree variant on the subgraphs
will prove to be the most practical ordering.
Anisotropy can have a significant e#ect on performance, both in terms of preconditioner
computing time and solution time. Our WND algorithm shows the most
promise for a high-performance algorithm that can exploit anisotropy, perhaps after
some tuning of the weights in the Laplacian matrix used. Robustness is still an issue;
we believe a more sophisticated weighting heuristic is necessary for further progress.


Appendix

. Testing data. Our test platform was a 180MHz Pentium Pro running
Windows/NT. We used MATLAB 5.1, with the algorithm for AINV written as a
MEX extension in C. Our AINV algorithm was a left-looking, column-by-column ver-
sion, with o#-diagonal entries dropped when their magnitude is below a user-supplied
tolerance and with the entries of D shifted to 10 -3 max |A| when their computed
magnitude We also make crucial use of the elimination tree;
in making a column conjugate with the previous columns, we only consider its descendants
in the elimination tree (the only columns that could possibly contribute
anything). This accelerates AINV considerably for low I F fill orderings-e.g., SHER-
MAN3 with AMDBAR ordering and a drop tolerance of 0.1 is accelerated by a factor
of four! An upcoming paper [7] will explore this more thoroughly.
To compare the orderings we selected several matrices, mostly from the Harwell-Boeing
collection. First we found the amount of true I F fill caused by each ordering,
given in Table 4. We then determined drop-tolerances for AINV to produce preconditioners
with approximately N and 2N nonzeros, where N is the number of nonzeros
in the given matrix. For each matrix, ordering, and fill level we attempted to solve
using BiCGStab (CG for SPD matrices), where b was chosen so that the
correct x is the vector of all 1s. Tables 5, 6, and 7 give the CPU time taken for
preconditioner computation, the iterations required to reduce the residual norm by
a factor of 10 -9 , and the CPU time taken by the iterations. We halted after 1800
iterations; the daggers in Tables 6 and 7 indicate no convergence at that point.
In each box of the tables, the upper line corresponds to the unweighted ordering
and the lower line its weighted counterpart. In Tables 5, 6, and 7 the numbers on the
left of the box correspond to the low-fill tests and those on the right to the high-fill
tests.
To highlight the winning ordering for each matrix, we have put the best numbers
in underlined boldface.
878 ROBERT BRIDSON AND WEI-PAI TANG

Table
Comparison of I F fill caused by di#erent orderings. Nonzero counts are given in thousands
of nonzeros. In each box the upper number corresponds to the unweighted ordering, and the lower
number corresponds to its weighted counterpart.
Given Red-black AMDBAR ND MIP
Name n NNZ OutIn OutIn/* OutIn/* WND OutIn/*
1486 1506 419 1526 709
4468 4435 1615 113777 1736
28974 19841 6815 26054 12750
128 91 37 114 59
SHERMAN5 3312 21 1340 1122 414 334 465
432 430 54 833 52
ORDERING, ANISOTROPY, AND APPROXIMATE INVERSES 879

Table
Comparison of CPU time for preconditioner computation. In each box the upper numbers
correspond to the unweighted ordering, and the lower numbers correspond to its weighted counterpart.
The numbers on the left refer to the low-fill test, and the numbers on the right refer to the high-fill
test.
Given Red-black AMDBAR ND MIP
Name OutIn OutIn/* OutIn/* WND OutIn/*
1.7 1.6 2.0 1.8 1.3 1.5 2.5 3.2 1.3 1.5
ADD32 33.9 45.6 32.7 43.5 3.4 3.6 3.2 3.4 3.7 3.4
7.1 8.0 6.1 6.3 3.5 3.4 2.9 3.0 3.7 3.4
ALE3D 14.3 26.6 13.0 24.1 9.4 16.1 6.4 11.0 15.9 28.3
BCSSTK14 6.2 11.1 6.0 10.7 2.3 3.7 2.0 3.3 3.0 4.7
5.6 8.4 5.6 8.5 2.2 3.6 3.6 5.6 3.5 5.7
70.6 76.6 81.0 84.5 61.6 61.9 54.7 55.8 65.1 58.2
NASA1824 4.6 7.8 4.0 7.1 1.3 2.1 1.5 2.3 1.8 2.8
3.9 6.0 3.9 6.1 1.4 2.2 1.9 3.0 1.7 2.8
8.4 15.2 8.5 14.7 3.0 5.0 2.4 3.8 3.5 5.7
8.4 14.5 8.5 14.7 3.1 5.1 2.8 4.5 4.2 6.8
1.3 1.7 1.1 1.8 0.8 1.1 1.0 1.4 0.9 1.3
ORSREG1 7.6 10.7 4.7 6.4 2.8 3.9 2.7 3.8 3.3 4.6
6.4 8.8 5.3 7.1 2.9 3.6 4.6 6.4 4.4 6.1
PORES2 2.8 4.0 2.5 3.8 1.3 1.9 0.7 1.1 1.3 1.8
1.7 2.5 2.0 2.9 1.0 1.5 1.2 1.8 1.6 2.4
SAYLR4 8.5 12.1 5.9 7.1 2.8 3.8 2.3 2.8 4.1 5.0
4.6 6.1 4.1 4.8 2.9 3.5 2.3 2.6 3.3 3.9
SHERMAN2 5.8 11.1 5.4 9.9 3.0 5.0 3.0 5.0 3.2 5.3
5.2 8.9 5.0 8.8 3.2 5.5 4.0 6.8 3.0 5.1
SHERMAN5 12.1 21.2 9.5 16.4 4.1 6.6 3.6 5.7 4.1 6.3
7.6 11.7 7.3 11.6 4.3 6.7 4.8 7.3 4.1 6.3
SWANG1 18.7 29.4 18.2 28.0 4.2 5.8 3.2 4.4 6.4 9.2
18.2 29.8 14.8 23.3 4.1 4.2 3.4 4.6 5.9 8.5
WANG1 19.6 31.3 10.8 16.4 6.2 8.8 6.2 8.9 9.2 13.2
16.2 24.2 12.1 18.0 6.4 9.2 6.2 8.8 9.2 13.3
2.1 2.9 2.0 2.8 0.7 0.8 1.0 1.2 0.8 0.9
7.2 11.3 3.8 6.0 2.0 3.0 1.8 2.5 2.7 4.0
4.8 7.7 3.8 5.8 2.2 3.2 1.8 2.5 2.4 3.6
7.4 12.0 4.2 6.6 2.1 3.0 1.8 2.6 2.8 3.9
5.1 7.5 3.7 5.5 2.7 3.9 1.9 2.6 2.5 3.

Table
Comparison of iterations required to reduce residual norm by 10 -9 . In each box the upper
numbers correspond to the unweighted ordering, and the lower numbers correspond to its weighted
counterpart. The numbers on the left refer to the low-fill test, and the numbers on the right refer to
the high-fill test.
Given Red-black AMDBAR ND MIP
Name OutIn OutIn/* OutIn/* WND OutIn/*
28 21 34
19 9
28 37 21 33
43 26 28 42 19
PORES3 87 28 90 23 81 31 36 21
37 26 36 23 36 22 35 23 36 22
28 19 28 21
28 19
43
28 6 28 4
ORDERING, ANISOTROPY, AND APPROXIMATE INVERSES 881

Table
Comparison of time taken for iterations. In each box the upper numbers correspond to the
unweighted ordering, and the lower numbers correspond to its weighted counterpart. The numbers
on the left refer to the low-fill test, and the numbers on the right refer to the high-fill test.
Given Red-black AMDBAR ND MIP
Name OutIn OutIn/* OutIn/* WND OutIn/*
ALE3D 18.4 11.5 10.2 9.3 5.5 8.6 10.0 8.7 4.9 5.0
4.9 5.8 3.6 5.4 4.2 4.2 5.0 6.9 6.1 7.7
9.3 9.0 9.2 8.9 8.6 7.5 13.0 15.4 7.9 7.1
13.4 16.0 13.2 17.0 8.3 7.7 10.3 10.5 7.7 7.5
14.7 11.4 10.4 6.8 22.0 5.7 15.0 4.3 11.1 4.5
NASA1824 71.4 73.5 68.2 73.5 63.3 52.1 53.2 50.4 65.3 52.0
NASA2146 10.3 10.2 11.6 12.8 7.8 6.3 7.9 5.6 7.2 6.9
11.6 7.4 10.9 7.4 9.2 6.7 9.4 6.0 7.9 6.1
ORSREG1 2.4 1.9 2.3 1.4 2.3 1.3 2.2 2.0 2.4 1.4
2.2 1.5 2.1 1.2 2.2 1.5 2.0 1.3 2.3 1.6
PORES3 1.1 0.5 1.2 0.4 1.1 0.5 0.5 0.3 0.4 0.3
4.3 4.6 4.0 3.8
4.8 5.0 4.0 4.4 4.0 3.9 4.3 3.7 4.1 3.8
SAYLR4 71.8 4.6 3.9 6.4 4.7 4.6 8.9 4.3 8.3 4.3
65.6 4.8 6.8 4.6 6.2 4.4 4.8 4.3 4.9 4.3
SHERMAN5 2.8 2.5 2.5 2.4 2.6 2.0 2.3 2.0 2.4 2.2
2.9 2.3 2.9 2.3 2.5 2.0 2.7 2.1 2.4 2.0
3.2 3.4 2.9 2.6 3.1 2.7 2.9 2.9 3.0 3.0
3.2 3.2 3.1 2.8 3.0 2.6 3.0 2.8 3.0 2.8
2.3 0.3 0.4 0.3 0.6 0.3 0.3 0.



--R

An approximate minimum degree ordering algorithm
A sparse approximate inverse preconditioner for the conjugate gradient method
A sparse approximate inverse preconditioner for nonsymmetric linear systems
Numerical experiments with two approximate inverse preconditioners
Orderings for factorized sparse approximate inverse preconditioners

Refined algorithms for a factored approximate inverse.
Approximate inverse techniques for general sparse matrices
Approximate inverse preconditioners via sparse-sparse iterations
Spectral ordering techniques for incomplete LU preconditioners for CG methods
Weighted graph based ordering techniques for preconditioned conjugate gradient methods
Towards a cost-e#ective ILU preconditioner with high level fill
Decay rates for inverses of band matrices

An algebraic approach to connectivity of graphs
Improving the performance of parallel factorized sparse approximate inverse precon- ditioners
Computer Solution of Large Sparse Positive Definite Systems
The evolution of the minimum degree ordering algorithm
Predicting structure in sparse matrix computations
Parallel preconditioning with sparse approximate inverses
A fast and high quality multilevel scheme for partitioning irregular graphs
Factorized sparse approximate inverse preconditionings I.
The role of elimination trees in sparse factorization
Toward an e
--TR

--CTR
E. Flrez , M. D. Garca , L. Gonzlez , G. Montero, The effect of orderings on sparse approximate inverse preconditioners for non-symmetric problems, Advances in Engineering Software, v.33 n.7-10, p.611-619, 29 November 2002
Michele Benzi, Preconditioning techniques for large linear systems: a survey, Journal of Computational Physics, v.182 n.2, p.418-477, November 2002
preconditioner;conjugate gradient-type methods;anisotropy;approximate inverse;ordering methods