A Priori Sparsity Patterns for Parallel Sparse Approximate Inverse Preconditioners.
Parallel algorithms for computing sparse approximations to the inverse of a sparse matrix either use a prescribed sparsity pattern for the approximate inverse or attempt to generate a good pattern as part of the algorithm.  This paper demonstrates that, for PDE problems, the patterns of powers of sparsified matrices (PSMs) can be used a priori as effective approximate inverse patterns, and that the additional effort of  adaptive sparsity pattern calculations may not be required.  PSM patterns are related to various other approximate inverse sparsity patterns through matrix graph theory and heuristics concerning the PDE's Green's function.  A parallel implementation shows that PSM-patterned approximate inverses are significantly faster to construct than approximate inverses constructed adaptively, while often giving preconditioners of comparable quality.
Introduction
. A sparse approximate inverse approximates the inverse of a
(usually sparse) matrix A by a sparse matrix M . This can be accomplished, for
example in the least-squares method, by minimizing the matrix residual norm 1
F
with the constraint that M is sparse. In general, the degrees of freedom of this
problem are the nonzero values in M as well as their locations. A minimization
that considers all these variables simultaneously, however, is very complex, and thus
a simple approach is to prescribe the set of nonzeros, or the sparsity pattern S, of
before performing the minimization. The objective function (1.1) can then be
decoupled as the sum of squares of the 2-norms of the n individual columns, i.e.,
in which e j and m j are the jth columns of the identity matrix and of the matrix M ,
respectively. Each least-squares matrix is small, having a number of columns equal
to the number of nonzeros in its corresponding m j . If A is nonsingular, then the
least-squares matrices have full rank.
Thus the approximate inverse can be constructed by solving n least-squares problems
in parallel. However, sparse approximate inverses are attractive for parallel pre-conditioning
primarily because the preconditioning operation is a sparse matrix by
vector product. The cost of constructing the approximate inverse for a large matrix
is usually so high, especially with the adaptive pattern selection strategies described
below, that they are only competitive if they are constructed in parallel.
For diagonally dominant A, the entries in A 1 decay rapidly away from the diagonal
[17], and a banded pattern for M will produce a good approximate inverse.
This work was performed under the auspices of the U.S. Department of Energy by the Lawrence
Livermore National Laboratory under Contract W-7405-ENG-48.
y Center for Applied Scientic Computing, Lawrence Livermore National Laboratory, L-560, Box
808, Livermore, CA 94551 (echow@llnl.gov).
1 The form for the right approximate inverse is used here, which is notationally slightly clearer. If
the matrix is distributed over parallel processors by rows, a left approximate inverse can be computed
row-wise.
In a general setting, without application-specic information, it is not clear how best
to choose a sparsity pattern for M . Algorithms have been developed that rst compute
an approximate inverse with an initial pattern S; then S is updated and a new
minimization problem is solved either exactly or inexactly. This process is repeated
until a threshold on the residual norm has been satised, or a maximum number of
nonzeros has been reached [12, 15, 22]. We refer to these as adaptive procedures.
One such procedure is to use an iterative method, such as minimal residual,
starting with a sparse initial guess [12] to approximately minimize (1.2), i.e., nd
sparse approximate solutions to
Suppose a sparse initial guess for m j is used. The rst few iterates will be sparse.
To maintain sparsity, a strategy to drop small elements is usually used, either in
the search direction or the iterates. No prescribed pattern S is necessary since the
sparsity pattern emerges automatically. For this method to be e-cient, sparse-sparse
operations must be used: the product of a sparse matrix by a sparse vector with p
nonzeros only involves p columns of the sparse matrix.
Another adaptive procedure, called SPAI [15, 22], uses a numerical test to determine
which nonzero locations should be added to the current sparsity pattern. For
the jth column, the numerical test for adding a nonzero in location k has the form
tolerance
is the residual for a given sparsity pattern of column j, and m j is
the current approximation. The test (1.4) is a lower bound on the improvement in the
square of the residual norm when the pattern of m j is augmented. Entry k is added
if the tolerance is satised, or if the left-hand side of (1.4) is large compared to its
values for other k. The cost of performing this test is a sparse dot product between r
and column k of A for each location to test. Interprocessor communication is needed
to test a k corresponding to a column not on a local processor.
These adaptive algorithms utilize the additional degrees of freedom in minimizing
aorded by the locations of the nonzeros in M and have allowed much more
general problems to be solved than before. Adaptive methods, however, tend to be
very expensive. Thus, this paper focuses on the problem of selecting S in a preprocessing
step so that a sparse approximate inverse can be computed immediately by
minimizing (1.1). Section 2 rst examines the sparsity patterns that are produced by
both non-adaptive and adaptive schemes. Section 3 tests the idea of using the patterns
of powers of sparsied matrices (PSM's) as a priori sparsity patterns for sparse
approximate inverses. Numerical tests are presented in Section 4, with comparisons
to both sequential and parallel versions of some current methods. Finally, Section 5
draws some conclusions.
2. Graph interpretations of approximate inverse sparsity patterns.
2.1. Use of the pattern of A and variants. The structure of a sparse matrix
A of order n is the directed graph G(A), whose vertices are the integers
whose edges (i ! j) correspond to nonzero o-diagonal entries in A. (This notation
usually implies matrices with all nonzero diagonal entries.) A subset of G(A) is a
directed graph with the same vertices, but with a subset of the edges in G(A). The
A PRIORI PATTERNS FOR SPARSE APPROXIMATE INVERSES 3
graph G(A) is a representation of the sparsity pattern of A, and when it is clear from
the context, we will not distinguish between them.
The structure of a vector x of order n is the subset of ng that corresponds
to the nonzero entries in x. When there is an associated matrix of order n, we often
refer to the structure of x as a subset of the vertices of the associated matrix. Notice
then, that the structure of column j of a matrix A is the set of vertices in G(A) that
have edges pointing to vertex j plus vertex j itself. The structure of row j is vertex
plus the set of vertices pointed to by vertex j.
The inverse of a matrix shows how each unknown in a linear system depends
on the other unknowns. The structure of the matrix A shows only the immediate
dependencies. This suggests that in the structure of A 1 there is an edge (i !
whenever there is a directed path from vertex i to vertex j in G(A) [21] (if A
is nonsingular, and ignoring coincidental cancellation). This structure is called the
transitive closure of G(A), and is denoted G  (A). For an irreducible matrix, this result
says that the inverse is a full matrix, but it does suggest the possibility of truncating
the transitive closure process to approximate the inverse by a sparse matrix.
A heuristic that is often employed is that vertices closer to vertex j along directed
paths are more important, and should be retained in an approximate inverse sparsity
pattern. This idea is supported by the decay in the elements observed by Tang [30]
in the discrete Green's function for many problems.
These sparsity patterns were rst used by Benson and Frederickson [4] in the
symmetric case, who also dened matrices with these patterns to be q-local matrices.
Given a graph G(A) of a structurally symmetric matrix A with a full diagonal, the
structure of the jth column of a q-local matrix consists of vertex j and its qth level
nearest-neighbors in G(A). A 0-local matrix is a diagonal matrix, while a 1-local
matrix has the same sparsity pattern as A.
The sparsity pattern of A is the most common a priori pattern used to approximate
A 1 . It gives good results for many problems, but can usually be improved, or
fails for many other problems. One improvement is to use higher levels of q. Unfortu-
nately, the storage for these preconditioners grows very quickly when q is increased,
and even impractical in many cases [20].
Huckle [24] proposed similar patterns which may be more eective when A is
nonsymmetric. These include the patterns corresponding to the graphs G((A T
The density of the former, in
particular, grows very quickly with increasing k. Primarily, these patterns are useful
as envelope patterns from which the adaptive SPAI algorithm can select its pattern.
This gives an upper bound on the interprocessor communication required by a parallel
implementation [23].
Cosgrove and Daz [14] proposed augmenting the pattern of A without going
to the full 2-local matrix. They suggested adding nonzeros to m j in a way that
minimizes the number of new rows introduced into the jth least-squares matrix (in
expression (1.2)). The augmented structure is determined only from the structure
of A. Kolotilina and Yeremin [28] proposed similar heuristics for augmenting the
sparsity pattern for factorized sparse approximate inverses.
Sparsication. Instead of augmenting the pattern of A, it is also possible to
diminish the pattern of A when A is relatively full. This can be accomplished by spar-
sication (dropping small elements in the matrix A) and using the resulting pattern.
This was introduced by Kolotilina [26] for computing sparse approximate inverses for
dense matrices (see also [10, 32]), and Kaporin [25] for sparse matrices and factorized
approximate inverses.
Sparsication can be combined with the use of higher level neighbors. Tang [30]
showed that sparsifying a matrix prior to applying the adaptive SPAI algorithm is
eective for anisotropic problems. The observation is that the storage and therefore
operation count required for preconditioners produced this way are much smaller.
This technique can generate patterns that are those of powers of sparsied matrices.
The idea of explicitly combining sparsication with the use of higher level neighbors
was used by Alleon et al. [1], who attributes the technique to Cosnuau [16]. For
approximating the inverse of dense matrices in electromagnetics, however, their tests
showed that higher levels were not warranted. Tang and Wan [31] also used a sparsi-
cation before applying a q-local matrix pattern, for q > 1, for approximate inverses
used as multigrid smoothers. They showed that the sparsication does not cause a
deterioration in convergence rate for their problems. Both the work by Alleon et al.
and Tang and Wan represent the rst uses of PSM patterns.
Instead of applying the sparsication to A, it is also appropriate in some cases to
apply the sparsication to the sparse approximate inverse after it has been computed
[27, 31]. This is useful to reduce the cost of using the approximate inverse when it is
relatively full.
2.2. Insights from adaptive schemes. Adaptive schemes can generate patterns
that are very dierent from the pattern of A, for example, the generated patterns
can be much sparser than A. Nevertheless, the patterns produced by adaptive schemes
can be interpreted using the graph of A.
Consider rst the approximate minimization method described in Section 1. The
following algorithm nds a sparse approximate solution to
minimal residual iterations. A dropping strategy for elements in the search direction
r is encapsulated by the function \drop" in step 4, which may depend on the current
pattern of m.
Algorithm 2.1. Sparse approximate solution to
1. m := sparse initial guess
2. r := e j Am
3. Loop until krk < tol or reached max. iterations
4. d := drop(r)
5. q := Ad
6.  := (r;q)
7. m := m+ d
8. r := r q
9. EndLoop
The elements in r not already in m are candidates for new elements in m. The
vector r is generated essentially by the product Am and thus the structure of r is the
set of vertices that have edges pointing to vertices in the structure of m. If the initial
guess for m consists of a single nonzero element at location j, then the structure of
grows outward from vertex j in G(A) with each iteration of the above algorithm.
If the search direction in the iterative method is A T r instead of r, then the kth
entry in the search direction is r T Ae k . A dropping strategy based on the size of the
these entries is similar to one based on the test (1.4) which attempts to minimize the
updated residual norm. In this method, the candidates for new elements are the rst
and second level neighbors of the vertices in m (for nonsymmetric A, the directions
of the edges are important).
A PRIORI PATTERNS FOR SPARSE APPROXIMATE INVERSES 5
This is exactly the graph interpretation for the SPAI algorithm (see Huckle [24]).
When computing column j of M , vertices far from vertex j will not enter into the
pattern, at least initially. Algebraically, this means that the nonzero locations of r
and Ae k do not intersect, and the value of the test is zero. An e-cient implementation
of SPAI uses these graph ideas to narrow down the indices k that need to be checked.
An early parallel implementation of SPAI [18] tested only the rst level neighbors
of a vertex, rather than both the rst and second levels. This is a good approximation
in many cases. This implementation also assumed A is structurally symmetric, so
that one-sided interprocessor communication is not necessary. A more recent parallel
implementation of SPAI [2, 3] implements the algorithm exactly. This code implements
one-sided communication with the Message Passing Interface (MPI), and uses
dynamic load balancing in case some processors nish computing their rows earlier
than others.
3. Patterns of powers of sparsied matrices.
3.1. Graph interpretation. In Section 2, we observed that prescribed sparsity
patterns or patterns generated by the adaptive methods are generally subsets of the
pattern of low powers of A (given that A has a full diagonal) and typically increase in
accuracy with higher powers. Clearly all these methods are related to the Neumann
series or characteristic polynomial for A [24].
The structure of column j of the approximate inverse of A is a subset of the
vertices in the level sets (with directed edges) about vertex j in G(A). Good vertices
to choose are those in the level sets in the neighborhood of vertex j, but the algorithms
dier in how these vertices are selected.
For convection-dominated and anisotropic problems, upstream vertices or vertices
in the preferred directions will have a greater in
uence than others on column j of
the inverse. Figure 3.1 shows the discrete Green's function for a point on a PDE
with convection. The nonsymmetry of the function shows that upstream nonzeros in
a row or column of the exact inverse are greater in magnitude than others. Without
additional physical information such as the direction of
ow, however, it is often
possible to use sparsication to identify the preferred or upstream directions.5152501020300.10.30.50.7
Fig. 3.1. Green's function for a point on a PDE with convection.
We examined the sparsity patterns produced by the adaptive algorithms and tried
to determine if they could have been generated by simpler graph algorithms. For some
simple examples, it turned out that the structures produced are exactly or very close
to the transitive closures of a subset of G(A), i.e., of the structure of a sparsied
6 EDMOND CHOW
matrix. In Figure 3.2 we show the structures of several matrices: (a) ORSIRR2 from
the Harwell-Boeing collection, (b) A 0 , a sparsication of the original matrix, (c) the
transitive closure G  (d) the structure produced by the SPAI algorithm. This
latter gure was selected from [2] which shows it as an example of an eective sparse
approximate inverse pattern for this problem. (There are, however, some bothersome
features of this example: the approximate inverse is four independent diagonal blocks.)
Note that we can approximate the adaptively generated pattern (d) very well by the
pattern (c) generated using the transitive closure.
(a)
(b) Sparsied ORSIRR2
(c) Transitive closure of (b)
(d) Pattern from SPAI
Fig. 3.2. The adaptively generated pattern (d) can be approximated by the transitive closure (c).
3.2. How to sparsify. The simplest method to sparsify a matrix is to retain
only those entries in a matrix greater than a global threshold, thresh. In the example
of

Figure

was used. It was important, however, to make sure that
the diagonal elements were retained, otherwise a structurally singular matrix would
have resulted in this case. In general, the diagonal should always be retained.
A PRIORI PATTERNS FOR SPARSE APPROXIMATE INVERSES 7
One strategy for choosing a threshold is to choose one that retains, for example,
one-third of the original nonzeros in a matrix. Fewer nonzeros should be retained
if powers of this sparsied matrix have numbers of nonzeros that grow too quickly.
This is how thresholds were chosen for the small problems tested in Section 4. The
number of levels used may be increased until a preconditioner reaches a target number
of nonzeros. The best choices for these parameters will be problem-dependent. For
special problems, this strategy may not be eective, for example, when a matrix
contains only a few unique values.
When a matrix is to be sparsied using a global threshold, how the matrix is
scaled becomes important. It is often the case that a matrix contains many dierent
types of equations and variables that are not scaled the same way. For example,
consider the matrix@
which has its rst row and column scaled by a large number, Z. If a threshold Z is
chosen, and if the diagonal of the matrix is retained, the third row of the sparsied
matrix has become independent of the other rows. We thus apply the thresholding to
a matrix that has been symmetrically scaled so that it has all ones on its diagonal,
e.g., for the above matrix, the scaled matrix is@
A threshold less than or equal to 1 will guarantee that the diagonal is retained. The
scaling also makes it easier to choose a threshold. This method of scaling is not
foolproof, but does avoid some simple problems.
In the graph of the sparsied matrix A 0 , each vertex should have some connections
to other vertices. This can be accomplished by sparsifying the matrix A such that
one retains at least a xed number of edges to or from each vertex, for example,
the ones corresponding to the largest matrix values. Given a parameter ll, this can
be implemented for column j by selecting the diagonal (jth) element plus the ll 1
largest o-diagonal elements in column j of the original matrix. This guarantees that
there are ll 1 vertices with edges into vertex j. Applied row-wise, this guarantees
ll 1 vertices with edges emanating from vertex j. Again, we choose explicitly to
keep the diagonal of the matrix; thus each column (or row) has at least ll nonzeros.
Choosing ll may be simpler and more meaningful than choosing a threshold on the
matrix values. Dierent values of ll may be used for dierent vertices, depending on
the vertex's initial degree (number of incident edges).
denote the structure of a matrix A 0 that has been sparsied from A.
denote the structure (with the same set of vertices) that has an edge
whenever there is a path of distance i or less in S 0 . The structure S i 1 is a subset
of S i . In matrix form, S
ignoring coincidental cancellation. These are
called level set expansions of a sparsied matrix or patterns of powers of a sparsied
matrix (PSM patterns). Heuristic 3 tested by Alleon et al. [1] is equivalent to S i using
a variable ll at each level to perform sparsication.
We mention that it is also possible to perform sparsication on S i after every
level set expansion. We denote this variant by S
i . For this variant, values need to be
computed, and we propose the following, which stresses the larger elements in A. If
\drop" denotes a sparsication process, then we can dene A
and S
note that S
1 is not generally a subset of
. More complicated strategies are possible; the thresholds can be dierent for each
level i. Note that determining S
i is much more di-cult than determining S i since
values need to be computed.
3.3. Factorized forms of the approximate inverse. Sparse approximate
inverses for the Cholesky or LU factors of A are often used. The analogue of the
least-squares method (minimization of (1.1)) here is the factorized sparse approximate
inverse (FSAI) technique of Kolotilina and Yeremin [28], implemented in parallel by
Field [19]. If the normal equations method is used to solve the least-squares systems,
the Cholesky or LU factors are not required to compute the approximate inverse. This
means, however, that the adaptive pattern selection schemes cannot be used, since
the matrix whose inverse is being approximated is not available. A priori sparsity
patterns must be used instead.
Given sparse approximate inverse approximates U 1 and
sparse matrices G and H , respectively, so that
The patterns for G and H should be chosen such that the pattern of GH is close
in some sense to good patterns for approximating A 1 . Supposing that S is a good
pattern for A 1 , then the upper and lower triangular parts of S are good patterns
for G and H , since the pattern of GH includes the pattern S. These patterns will be
tested in Section 4.
It may also be possible to use the patterns of the powers of the exact or approximate
L and U if they are known. These L and U factors are not discretizations of
PDE's, but their inverses are often banded with elements decaying rapidly away from
the main diagonal. This technique may be appropriate if approximate L and U are
available, for example from a very sparse incomplete LU factorization.
As opposed to the inverses of irreducible matrices, the inverses of Cholesky or
LU factors are often sparse. An ordering should thus be applied to A that gives
factors whose inverses can be well approximated by sparse matrices. Experimentally,
fewer nonzeros in the exact inverse factors translates into lower construction cost and
better performance for factorized approximate inverses computed by an incomplete
biconjugation process [7, 8]. The transitive closure can be used to compute the number
of nonzeros in the exact inverse of a Cholesky factor, based on the height of all
the nodes in the elimination tree. This has lead to reordering strategies that approximately
minimize the height of the elimination tree and thus the number of nonzeros in
the inverse factors, and allows some prediction of how well these approximate inverses
might perform on a given problem [7, 8].
3.4. Approximate inverse of a Schur complement. To determine a good
pattern for a Schur complement matrix, we notice that
is the Schur complement. Thus a good sparsity pattern for
S 1 can be determined from a good sparsity pattern for A 1 ; it is simply the (2,2)
A PRIORI PATTERNS FOR SPARSE APPROXIMATE INVERSES 9
block of a good sparsity pattern for A 1 . B should be of small order compared to the
global matrix or else the method will be overly costly. In a code, it may be possible
to compute the approximate inverse of A and extract the approximation to S 1 , or
compute a partial approximate inverse, i.e., those rows or columns of the approximate
inverse that correspond to S 1 [11]. Again, S should be of almost the same order as
A.
3.5. Parallel computation. Computation of S i is equivalent to structural sparse
matrix-matrix products of sparsied matrices. The computation can also be viewed as
n level set expansions, one for each row or column, which can be performed in parallel.
For vertices that are near other vertices on a dierent processor, some communication
will be necessary. Communication can be reduced by partitioning the graph of the
sparsied matrix among the processors such that the number of edge-cuts is reduced.
Unfortunately, in general, one-sided communication is required to compute sparse
approximate inverses. Processors need to request rows from other processors, and a
processor cannot predict which rows it will need to send. One-sided communication
may be implemented in MPI by having each processor occasionally probe for messages
from other processors. The latency between probes is a critical performance factor
here. In a multithreaded environment, it is possible to dedicate some threads on
a local processing node to servicing requests for rows (server threads), while the
remaining threads compute each row and make requests for rows when necessary
(worker threads).
Consider a matrix A and an approximate inverse M to be computed that are
partitioned the same way by rows across several processors. Algorithm 3.1 describes
one organization of the parallel computation. Each processor computes a level set
expansion for all of its rows before before continuing on to the next level. At each
level, the requests and replies to and from a processor are coalesced, allowing fewer
and larger messages to be used. Like in [3], external rows of A are cached on a
processor in case they are needed to compute other rows. There is no communication
during the numerical phase when the values of M are being computed. This algorithm
was implemented using occasional probing for one-sided communication.
Algorithm 3.1. Parallel level set expansions for computing S i
Communicate rows
1. Initialize the set of vertices V to empty
2. Sparsify all the rows on the local processor
3. Merge the structures of all the locally sparsied rows into V
4. For level
5. For nonlocal k 2 V, request and receive row k
6. Sparsify received rows
7. Merge structures of new sparsied rows into V
8. EndFor
9. For nonlocal k 2 V, request and receive row k
Compute structure of each row
10. For each local row j
11. Initialize V j to a single entry in location j
12. For level
13. For new sparsied structure of row k into V j
14. EndFor
15. EndFor
Compute values of each row
16. For each local row j, nd
has the pattern V j
We also implemented a second parallel code, which has the following features:
multithreaded, to take advantage of multiple processors per shared-memory
node on symmetric multiprocessor computers
uses server and worker threads to more easily implement one-sided commu-
nication
uses a simpler algorithm than Algorithm 3.1: computes each row and performs
all the associated communications before continuing on to the next row; when
multiple threads are used, this avoids worker threads needing to synchronize
and coordinate which rows to request from other nodes; smaller messages are
used, but communication is also spread over the entire execution time of the
algorithm
scalable with its use of memory, but is thus also slower than the rst version
which used direct-address tables (traded memory for faster computation)
Timings for this second code will be reported in the next section. Some limited
timings for the rst code will also be shown. We are also working on a factorized
implementation for symmetric matrices, which will guarantee that the preconditioner
is also symmetric. This implementation makes a simple change in step 16 of Algorithm
3.1, and does not require one-sided communication when the full matrix is stored.
4. Numerical tests.
4.1. Preconditioning quality. First we test the quality of sparsity patterns
generated by powers of sparsied matrices on small problems from the Harwell-Boeing
collection. In particular, we chose problems that were tested with SPAI [22] in order
to make comparisons. We performed tests in exactly the same conditions: we solve the
same linear systems using GMRES(20) to a relative residual tolerance of 10 8 with a
zero initial guess. We report the number of GMRES steps needed for convergence, or
indicate no convergence using the symbol y. Right preconditioning was used.
In

Tables

4.1 to 4.5, we show test results for S i for both unfactored (column 2)
and factored (column 3) forms of the approximate inverse. We compare the results to
the least-squares (LS) method using the pattern of the original matrix A, and FSAI,
the least-squares method for the (nonsymmetric) factored form [28], again using the
pattern of the original matrix A. For the unfactored form, we also display the result
of the SPAI method reported in [22], using their choice of parameters. Adaptive
methods for factored forms are also available [5, 6, 29] but were not tested here.
Global thresholds (shown for each table) on a scaled matrix were used to perform
these sparsications. In the tables, we also show the number of nonzeros nnz in the
unfactored preconditioners (the entry for LS/FSAI is the number of nonzeros in A).
The results show that preconditioners of almost the same quality as the adaptive
SPAI can be achieved using the S i patterns. In some cases, even better preconditioners
can result, sometimes with even less storage (Table 4.1). The results also show that
using the pattern of A does not generally give as good a preconditioner for these
problems.
SHERMAN2 is a relatively hard problem for sparse approximate inverses. The
result reported in [22] shows that SPAI could reduce the residual norm by 10 5 (the
preconditioner with 26327 nonzeros in 7 steps. The results
are similar with S i patterns, but the full residual norm reduction can be achieved
with an approximate inverse that is denser. Note that in this case, the sparsication
A PRIORI PATTERNS FOR SPARSE APPROXIMATE INVERSES 11

Table
Iteration counts for ORSIRR2,
Pattern unfactored factored nnz
LS/FSAI 335 383 5970
SPAI 84 5318

Table
Iteration counts for SHERMAN1,
Pattern unfactored factored nnz
LS/FSAI 145 456 3750
SPAI
threshold was applied to the original matrix rather than to the diagonally scaled
matrix, although the matrix has values over 27 orders of magnitude; the diagonal
scaling is not foolproof. Factorized approximate inverses were not eective for this
problem with these patterns.
SAYLR4 is a relatively hard problem for GMRES. Grote and Huckle [22] report
that SPAI could not solve the problem with GMRES, but could with BiCGSTAB.
This is also true for S i patterns; the results in Table 4.5 are with BiCGSTAB.
There are, of course, many problems that are di-cult for PSM-patterned approximate
inverses. These include NNC666 and GRE1107 from the Harwell-Boeing
collection, and FIDAP problems from Navier-Stokes simulations. These problems,
however, can be solved using adaptive methods [12]. These problems pose di-culties
for PSM patterns because the Green's function heuristic is invalid; the problems either
are not PDE problems, or have been modied (e.g., the FIDAP problems used a
penalty formulation).
In

Tables

4.6 and 4.7, we show test results for S
i for the unfactored form of the
approximate inverse. S
0 and the LS patterns were the same for these problems. Since
is not generally a superset of S
there is no guarantee that S
is a better pattern
than S
in terms of the norm of R = I AM . To show this, we also display these
matrix residual norms. The parameter ll (shown for each table) was used to sparsify
these matrices after each level set expansion. Again, the results show that the a priori
methods can approach the quality of the adaptive methods very closely.
4.2. Parallel timing results. The results above show that the preconditioning
quality for PDE problems is not signicantly degraded by using the non-adaptive
schemes based on powers of sparsied matrices. In this section we illustrate the main
advantage of these preconditioners: their very low construction costs compared to the
adaptive schemes.
In this and Section 4.3 we show results using the multithreaded version of our
code, ParaSAILS (parallel sparse approximate inverse, least squares). Like the parallel
version of SPAI [3] with which we make comparisons, ParaSAILS is implemented

Table
Iteration counts for SHERMAN2,
Pattern unfactored factored nnz
LS/FSAI y y 23094
SPAI y 26327

Table
Iteration counts for PORES3,
Pattern unfactored factored nnz
LS/FSAI y y 3474
SPAI 599 16745
as a preconditioner object in the ISIS++ solver library [13]. Both these codes generate
a sparse approximate inverse partitioned across processors by rows; thus left
preconditioning is used. In all the codes, the least-squares problems that arose were
solved using LAPACK routines for QR decomposition. For problems with relatively
full approximate inverses, solving these least-squares problems takes the majority of
the computing time.
Tests were run on multiple nodes of an IBM RS/6000 SP supercomputer at the
Lawrence Livermore National Laboratory. Each node contains four 332 MHz PowerPC
CPU's. Timings were performed using user-space mode, which is much more
e-cient than internet-protocol mode. However, nonthreaded codes can only use one
processor per node in user-space mode. We tested SPAI with one processor per node,
and ParaSAILS with up to four processors per node. The iterative solver and matrix-vector
product codes were also nonthreaded, and used only one processor per node.
The rst problem we tested is a nite element model of three concentric spherical
shells with dierent material properties. The matrix has order 16881 and has
nonzeros. The SPAI algorithm using the default parameters (target
residual norm for each row < 0:4) produced a much sparser precondi-
tioner, with 171996 nonzeros, and solved the problem using GMRES(50) to a tolerance
steps. For comparison purposes, we chose parameters for ParaSAILS
that gave a similar number of nonzeros in the preconditioner. In particular, S 3 with
an ll parameter of 3 gave a preconditioner with 179550 nonzeros, and solved the
problem in 331 steps. Figure 4.1 shows the two resulting sparsity patterns. Table
4.8 reports, for various numbers of nodes (npes), the wall-clock times for the preconditioner
setup phase (Precon), the iterative solve phase (Solve), and the total time
(Total). The time for constructing the preconditioner in each code includes the time
for determining the sparsity pattern. Due to the relatively small size of this (and
the next) problem, only one worker and one server thread was used per node (i.e.,
two processors per node) in the ParaSAILS runs; one processor was used in the SPAI
runs.
For comparison, in Table 4.9 we show the results using the rst (nonthreaded,
occasional MPI probes for one-sided communication) version of the code. This code
A PRIORI PATTERNS FOR SPARSE APPROXIMATE INVERSES 13

Table
Iteration counts for SAYLR4,
Pattern unfactored factored nnz
LS/FSAI y y 22316
SPAI 67 84800

Table
Iteration counts for SHERMAN3,
steps nnz kRkF
LS y 20033 17.3620
SPAI 264 48480
is much faster because it uses direct-address arrays to quickly merge the sparsity
patterns of rows; direct-address arrays have length the global size of the matrix and
are not scalable.
(a) ParaSAILS (b) SPAI
Fig. 4.1. Structure of the sparse approximate inverses for the concentric shells problem.
We tested a second, larger problem which models the work-hardening of metal by
squeezing it to make it pancake-like In this particular
example, the pattern of a sparsied matrix (S 0 ) with an ll parameter of 3 lead to a
good preconditioner. ParaSAILS produced a preconditioner with 141848 nonzeros and
solved the problem in 142 steps; SPAI produced a preconditioner with 120192 nonzeros
14 EDMOND CHOW

Table
Iteration counts for SHERMAN4,
steps nnz kRkF
LS 199 3786 6.2503
SPAI 86 9276

Table
Timings for concentric shells problem.
ParaSAILS SPAI
npes Precon Solve Total Precon Solve Total
and solved the problem in 139 steps. Results for varying numbers of processing nodes
are shown in Table 4.10. The results show that the non-adaptive ParaSAILS algorithm
implemented here is many times faster than the adaptive SPAI algorithm.
4.3. Implementation scalability. In this section, we experimentally investigate
the implementation scalability of constructing sparse approximate inverses with
ParaSAILS. Let T (n; p) be the time to construct an approximate inverse of order n
on a parallel computer using p processors. We dene the scaled e-ciency to be
E(n; p)  T (n; 1)=T (pn; p):
If E(n; then the implementation is perfectly scalable, i.e., one
could double the size of the problem and the number of processors without increasing
the execution time. However, as long as E(n; p) is bounded away from zero for a xed
n as p is increased, we say that the implementation is scalable.
We consider the 3-D constant coe-cient PDE
in
@
discretized using standard nite dierences on a uniform mesh, with the anisotropic
parameters This problem has been used to test the
scalability of multigrid solvers [9]. The problems are a constant size per compute node,
from 10 3 to local problem sizes. Node topologies of 1 3 to 5 3 were used. Thus the
largest problem was over a cube with (60  unknowns. Each node
used all four processors (4 worker threads, 1 server thread) for this problem in the
preconditioner construction phase (i.e., 500 processors in the largest conguration).
A threshold for ParaSAILS was chosen so that only the nonzeros along the
strongest (z) direction are retained, and the S 3 pattern was used. Although this is a
symmetric problem, the preconditioner is not symmetric, and we used GMRES(50)
as the iterative solver with a zero initial guess. The convergence tolerance was 10 6 .
A PRIORI PATTERNS FOR SPARSE APPROXIMATE INVERSES 15

Table
Concentric shells problem: Timings for preconditioner setup using direct addressing.
npes Precon

Table
Timings for work-hardening problem.
ParaSAILS SPAI
npes Precon Solve Total Precon Solve Total

Table

4.11 shows the results, including wall-clock times for constructing the preconditioner
and the iterative solve phase, the number of iterations required for convergence,
the average time for one iteration in the solve phase, and E p and E s , the scaled e-cien-
cies for constructing the preconditioner and for one step in the solve phase. Figures
4.2 and 4.3 graph the scaled e-ciencies E p and E s , respectively. The implementation
seems scalable for all values of p that may be encountered. For comparison, in Table
4.12, we show results for SPAI using a 40 3 local problem size. Larger problems led to
excessive preconditioner construction times. Again, one processor per node was used
for SPAI.
5. Conclusions. This paper demonstrates the eectiveness of patterns of powers
of sparsied matrices for sparse approximate inverses for PDE problems. As opposed
to many existing methods for prescribing sparsity patterns, PSM patterns use
both the values and structure of the original matrix, and very sparse patterns can be
produced. PSM patterns allow simpler direct methods of constructing sparse approximate
inverse preconditioners to be used, with comparable preconditioning quality to
adaptive methods, but with signicantly less computational cost. The numerical tests
show that the additional eort of adaptive sparsity pattern calculations is not always
required.

Acknowledgments

. The author is indebted to Wei-Pai Tang, who was one of
the rst to use sparsication for computing approximate inverse sparsity patterns.
John Gilbert was instrumental in directing attention to the transitive closure of a
matrix, and motivating the possibility of nding good patterns a priori. Michele
Benzi made helpful comments and directed the author to [24]. The author is also
grateful for the ongoing support of Robert Clay, Andrew
Esmond Ng, Ivan Otero, Yousef Saad, and Alan Williams, and nally for the cogent
comments of the anonymous referees.

Table
Timings, iteration counts, and e-ciencies for ParaSAILS.
local problem size
npes N Precon Solve Iter Solve/Iter Ep Es
local problem size
npes N Precon Solve Iter Solve/Iter Ep Es
local problem size
npes N Precon Solve Iter Solve/Iter Ep Es
local problem size
npes N Precon Solve Iter Solve/Iter Ep Es
50  50  50 local problem size
npes N Precon Solve Iter Solve/Iter Ep Es
local problem size
npes N Precon Solve Iter Solve/Iter Ep Es

Table
Timings, iteration counts, and e-ciencies for SPAI.
local problem size
npes N Precon Solve Iter Solve/Iter Ep Es
A PRIORI PATTERNS FOR SPARSE APPROXIMATE INVERSES 17
1200.10.30.50.70.9Number of nodes
Scaled
efficiency
Fig. 4.2. Implementation scalability of ParaSAILS preconditioner construction phase.
1200.10.30.50.70.9Number of nodes
Scaled
efficiency
Fig. 4.3. Implementation scalability of one step of iterative solution.



--R


An MPI implementation of the SPAI preconditioner on the T3E
A portable MPI implementation of the SPAI preconditioner in ISIS
Iterative solution of large sparse linear systems arising in certain multidimensional approximation problems



An ordering method for a factorized approximate inverse pre- conditioner
Semicoarsening multigrid on distributed memory machines
On a class of preconditioning methods for dense linear systems from boundary ele- ments
Approximate inverse techniques for block-partitioned matrices




Etude d'un pr
Decay rates for inverses of band matrices
Parallel implementation of a sparse approximate inverse preconditioner


Predicting structure in sparse matrix computations
Parallel preconditioning with sparse approximate inverses


A preconditioned conjugate gradient method for solving discrete analogs of
Explicit preconditioning of systems of linear algebraic equations with dense matrices
Factorized sparse approximate inverse preconditionings.
Factorized sparse approximate inverse precondition- ings I
Iterative Methods for Sparse Linear Systems
Towards an e
Sparse approximate inverse smoother for multi-grid
Preconditioning for boundary integral equations
--TR

--CTR
Kai Wang , Sang-Bae Kim , Jun Zhang , Kengo Nakajima , Hiroshi Okuda, Global and localized parallel preconditioning techniques for large scale solid Earth simulations, Future Generation Computer Systems, v.19 n.4, p.443-456, May
Robert D. Falgout , Jim E. Jones , Ulrike Meier Yang, Conceptual interfaces in hypre, Future Generation Computer Systems, v.22 n.1, p.239-251, January 2006
Robert D. Falgout , Jim E. Jones , Ulrike Meier Yang, Pursuing scalability for hypre's conceptual interfaces, ACM Transactions on Mathematical Software (TOMS), v.31 n.3, p.326-350, September 2005
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,
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
Dennis C. Smolarski, Diagonally-striped matrices and approximate inverse preconditioners, Journal of Computational and Applied Mathematics, v.186 n.2, p.416-431, 15 February 2006
Edmond Chow, Parallel Implementation and Practical Use of Sparse Approximate Inverse Preconditioners with a Priori Sparsity Patterns, International Journal of High Performance Computing Applications, v.15 n.1, p.56-74, February  2001
Oliver Brker , Marcus J. Grote, Sparse approximate inverse smoothers for geometric and algebraic multigrid, Applied Numerical Mathematics, v.41 n.1, p.61-80, April 2002
Luca Bergamaschi , Giorgio Pini , Flavio Sartoretto, Computational experience with sequential and parallel, preconditioned Jacobi--Davidson for large, sparse symmetric matrices, Journal of Computational Physics, v.188
Michele Benzi , Miroslav Tma, A parallel solver for large-scale Markov chains, Applied Numerical Mathematics, v.41 n.1, p.135-153, April 2002
Anwar Hussein , Ke Chen, Fast computational methods for locating fold points for the power flow equations, Journal of Computational and Applied Mathematics, v.164-165 n.1, p.419-430, 1 March 2004
P. K. Jimack, Domain decomposition preconditioning for parallel PDE software, Engineering computational technology, Civil-Comp press, Edinburgh, UK, 2002
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
preconditioned iterative methods;graph theory;sparse approximate inverses;parallel computing