Parallel Preconditioning with Sparse Approximate Inverses.
A parallel preconditioner is presented for the solution of general sparse linear systems of equations. A sparse approximate inverse is computed explicitly and then applied as a preconditioner to an iterative method. The computation of the preconditioner is inherently parallel, and its application only requires a matrix-vector product. The sparsity pattern of the approximate inverse is not imposed a priori but captured automatically. This keeps the amount of work and the number of nonzero entries in the preconditioner to a minimum. Rigorous bounds on the clustering of the eigenvalues and the singular values are derived for the preconditioned system, and the proximity of the approximate to the true inverse is estimated. An extensive set of test problems from scientific and industrial applications provides convincing evidence of the effectiveness of this approach.
Introduction
We consider the linear system of equations
The work of M. Grote was supported by an IBM fellowship. The work of T. Huckle
was supported by a research grant of the Deutsche Forschungsgemeinschaft.
y Scientific Computing and Computational Mathematics, Bldg. 460, Stanford University,
Stanford, CA 94305 (grote@sccm.stanford.edu).
z Institut f?r Angewandte Mathematik und Statistik, Universit?t W-urzburg, D-97074
W-urzburg, Germany (huckle@vax.rz.uni-wuerzburg.d400.de).
Here A is a large, sparse, and nonsymmetric matrix. Due to the size of A,
direct solvers become prohibitively expensive because of the amount of work
and storage required. As an alternative we consider iterative methods such as
gmres, bcg, bi-cgstab, and cg applied to the normal equations [6]. Given
the initial guess x 0 , these algorithms compute iteratively new approximations
x k to the true solution b. The iterate xm is accepted as a solution if the
residual r In general, the convergence
is not guaranteed or may be extremely slow. Hence, the original problem (1)
must be transformed into a more tractable form. To do so, we consider a
preconditioning matrix M , and apply the iterative solver either to the right
or to the left preconditioned system
Therefore, M should be chosen such that AM (or MA) is a good approximation
of the identity. As the ultimate goal is to reduce the total execution
time, both the computation of M and the matrix-vector product My should
be evaluated in parallel. Since the matrix-vector product must be performed
at each iteration, the number of nonzero entries in M should not exceed that
in A.
The most successful preconditioning methods in terms of reducing the
number of iterations, such as incomplete LU factorizations or SSOR, are
notoriously difficult to implement on a parallel architecture, especially for
unstructured matrices. Indeed, the application of the preconditioner in the iteration
phase requires the solution of triangular systems at each step, which is
difficult to parallelize because of the recursive nature of the computation (see
[2], section 4.4.4). Our aim is to find an inherently parallel preconditioner,
which retains the convergence properties of incomplete LU .
A natural way to achieve parallelism is to compute an approximate inverse
M of A, such that AM ' I in some sense. The evaluation of My is
then easy to parallelize, and will be cheap if M is sparse. Although such a
sparse approximate inverse does not always exist, it often occurs in applications
that most entries in A \Gamma1 are very small [5], [1]. For instance, if the
problem results from the discretization of a partial differential equation, it is
generally meaningful to look for a sparse approximate inverse. Polynomial
preconditioners with are inherently parallel, but do not lead to
as much improvement in the convergence as incomplete LU (see [2], section
3.
A different approach is to minimize kAM \Gamma Ik. Yet we cannot confine the
whole spectrum of AM to the vicinity of 1, or else we would be minimizing the
1-norm or the 2-norm. This would be too expensive, since
solution as long as kAM \Gamma Ik ? 1. To ensure fast convergence at a reasonable
cost, we need to cluster most eigenvalues and singular values about 1, but give
a few outliers additional freedom. We achieve this by minimizing kAM \Gamma Ik
in the Frobenius norm. Moreover, this choice naturally leads to inherent
parallelism, because the columns m k of M can be computed independently of
one another. Indeed, since
the solution of (3) separates into n independent least squares problems
Thus, we can solve (4) in parallel and obtain
an explicit approximate inverse M of A. If M is sparse, (4) reduces to n
small least squares problems, which can be solved very quickly [9], [13]. The
difficulty lies in determining a good sparsity structure of the approximate
inverse, or else the solution of (4) will not yield an effective preconditioner.
Therefore, we seek a method that captures the sparsity pattern of the main
entries of A \Gamma1 automatically, yet at a reasonable cost. We start with a given
sparsity pattern, such as diagonal, and augment M progressively until the 2-
norm of the residual is small enough or a maximal sparsity has been reached.
The key to the algorithm lies in the strategy used to determine the locations
of the entries in M . Our selection criterion is simple and cheap to compute,
yet it yields an effective preconditioner without generating excessive fill-in.
The extensive set of difficult test problems we consider at the end shows that
our algorithm produces a sparse and effective preconditioner.
The computation of approximate inverses, based on minimizing (4), has
been proposed by several authors. Yeremin et al compute a factorized sparse
approximate inverse [11], [12], [13], but only consider fixed sparsity patterns.
Simon and Grote solve (4) explicitly, but only allow for a banded sparsity
pattern in M [9], [10]. The approach of Cosgrove, Diaz, and Griewank [4]
is similar to ours, but differs in the criteria used for augmenting M . Chow
and Saad [3] use an iterative method to compute an approximate solution of
(4). Their method automatically generates new entries, to which they apply
a dropping strategy to remove the excessive fill-in appearing in M .
In section 2 we introduce the spai algorithm, which computes a sparse
approximate inverse of A. In section 3 we derive theoretical properties of the
spectrum of the preconditioned system (2). In sections 4 and 5 we present
a wide range of numerical experiments using test matrices from engineering
and scientific applications.
2 Computation of the Sparse Approximate Invers

We shall first show how to compute a sparse approximate inverse M for a
given sparsity structure. The matrix M is the solution of the minimization
problem (4). Since the columns of M are independent of one another, we only
need to present the algorithm for one of them, and we denote it by m k . Now
let J be the set of indices j such that m k (j) 6= 0 . We denote the reduced
vector of unknowns m k (J ) by -
I be the set of indices i such
that A(i; J ) is not identically zero. This enables us to eliminate all zero rows
in the submatrix A(:; J ) . We denote the resulting submatrix A(I; J ) by -
A .
Similarly, we define - e . If we now set
that solving (4) for m k is equivalent to solving
min
for -
. The n 1 \Theta n 2 least squares problem (5) is extremely small because A
and M are very sparse matrices. If A is nonsingular, the submatrix -
A must
have full rank. Thus, the QR-decomposition of -
A is
where R is a nonsingular upper triangular n 2 \Theta n 2 matrix. If we let -
the solution of (5) is
We solve (7) for each
This yields an
approximate inverse M , which minimizes kAM \Gamma Ik F for the given sparsity
structure.
Our aim is now to improve upon M by augmenting its sparsity structure to
obtain a more effective preconditioner. To do so, we shall reduce the current
error that is reduce kAm We
recall that m k is the optimal solution of the least squares problem (4), and
we denote its residual by
is exactly the k-th column of A \Gamma1 and cannot be improved upon.
We now assume that r 6= 0 and demonstrate how to augment the set of indices
J to reduce krk 2 . Since A and m k are sparse, most components of r are
zero, and we denote by L the remaining set of indices ' for which r(') 6= 0 .
Typically L is equal to I, since - r does not have exact zero entries in finite
precision. But if I does not contain k, it must be included in L since r(k) is
then equal to \Gamma1. To every ' 2 L corresponds an index set N ' , which consists
of the indices of the nonzero elements of A('; :) that are not in J yet. The
potential new candidates that might be added to J are contained in
~
We must now select new indices j that will lead to the most profitable reduction
in krk 2 . To do so in a cheap but effective way, we consider for each
J the one-dimensional minimization problem
The solution of (10) is
2: (11)
For each j we compute the 2-norm ae j of the new residual r
given by (11):
2: (12)
There is at least one index
J such that r T Ae j 6= 0, which will lead to a
smaller residual in (12). Otherwise
which would imply that r(L) is zero, since A(L; :) has full rank. We note
that J [ ~
J contains the column indices of all nonzero elements of A(L; :),
and that J " ~
We reduce ~
J to the set of the most profitable indices
j with smallest ae j and add it to J . We note that equation (9) was also used
in [4] to estimate the reduction in the residual, and can already be found in
[8].
Using the augmented set of indices J , we solve the sparse least squares
problem (4) again. This yields a better approximation m k of the k-th column
of A \Gamma1 . We repeat this process for each until the residual satisfies
a prescribed tolerance or a maximum amount of fill-in has been reached in
. The numerical study in section 4 shows that this iterative procedure
captures the main entries of A \Gamma1 extremely well.
Every time we augment the set of nonzero entries in m k , we solve the
least squares problem exactly. We shall now demonstrate how one can easily
update the QR decomposition and greatly reduce the amount of work. We
recall that J is the current set and that ~
J is the set of new indices that will
be added to m k . We denote by ~
I the new rows, which correspond to the
nonzero rows of A(:; J [ ~
J ) not contained in I yet, and by ~
2 the
number of indices in ~
I and ~
J . Thus, we need to replace in (5) the submatrix
A by the larger submatrix A(I [ ~
To solve (5), we use the known
QR decomposition of -
A to update the QR decomposition of the augmented
A A(I; ~
I ~
I ~
This requires only the computation of the QR decomposition of B 2 . We
note that A( ~
I already contains the indices of all nonzero
entries present in columns J . We let
R
is a ~
Using (16) we rewrite (15) as
I ~
'' I n 2
~
R
This procedure enables us to add new indices to J and solve the least squares
problem for the optimal solution, without recomputing the full QR decomposition
at each step. It generalizes the updating strategy proposed in [4] for a
single entry, by allowing for several new entries at a time.
If we do not stop the process, the algorithm will compute the k-th column
of A \Gamma1 . In practice, however, we stop the process once the prescribed tolerance
is met or a maximal amount of fill-in has been reached. We now present
the full spai algorithm, where spai stands for SParse Approximate Inverse:
The SPAI Algorithm:
For every column m k of M :
(a) Choose an initial sparsity J .
(b) Compute the row indices I of the corresponding nonzero entries and the
QR decomposition (6) of A(I; J ) . Then compute the solution m k of
the least squares problem (4), and its residual r given by (8).
While
(c) Set L equal to the set of indices ' for which r(') 6= 0.
(d) Set ~
J equal to the set of all new column indices of A that appear
in all L rows but not in J .
(e) For each j 2 ~
J solve the minimization problem (10).
(f) For each j 2 ~
J compute ae j given by (12), and delete from ~
J all
but the most profitable indices.
(g) Determine the new indices ~
I and update the QR decomposition
using (17). Then solve the new least squares problem, compute
the new residual
I and
J .
Remarks:
1. The initial sparsity structure of M is arbitrary and may be chosen
empty or diagonal, if no a priori information about the sparsity of A \Gamma1
is available. Yet, to solve a sequence of problems with similar sparsity
patterns but varying entries in A, a clever initial guess for the initial
sparsity is to choose the sparsity of the previously computed approximate
inverse. This greatly reduces the computational cost of M , since
the initial sparsity structure would be almost optimal. In all our numerical
examples the initial sparsity of M was chosen diagonal.
2. In addition to the stopping criterion on krk 2 , we constrain the loop to
a maximal number of iterations to limit the maximal fill-in per column
in M . This threshold was almost never reached and the total number
of nonzero entries in M is usually comparable to the amount in A.
3. In (f) we first reduce ~
J to the set of indices j such that ae j is less than
or equal to the mean value of all ae j . From the remaining indices we
keep at most s indices with smallest ae j . Here s should be a small
integer to avoid excessive fill-in, and we have set s equal to 5 in most
numerical calculations. This criterion is very cheap to compute and
removes useless indices effectively. It is parameter-free, since it uses a
dynamic mean value criterion and does not require a threshold input
by the user. A more sophisticated weighting could be applied to the
distribution of the ae j to control the rate at which M is being filled.
4. When we update the QR decomposition (17), we store the Householder
matrices resulting from the factorizations of the matrices B 2 separately,
and never construct Q explicitly.
5. The selection process in (c) may be restricted to the most easily accessible
rows to minimize communications and data flow. It can also be
restricted only to the largest elements in r.
6. The one-dimensional minimization can be replaced by another minimization
method such as steepest descent or the exact minimization
problem related to J [ fjg. From our experience, the former is not
accurate enough, and the latter is too expensive.
7. The approximate inverse M computed with the spai algorithm is permutation
invariant. If A is replaced by P 1 AP 2 , where P 1 and P 2 are
permutation matrices, we obtain P T
1 instead of M .
The spai algorithm may also be applied to compute a sparse approximate
left inverse of A, which can be used as a left preconditioner in (2). This may
yield a better result if all the rows of A \Gamma1 are sparse but a few columns in
A are full. Such a case is discussed in section 5.
If the iterative solver preconditioned with the current M does not con-
verge, it is easy to improve upon M using the sparsity of M as initial sparsity
for the spai algorithm. The iteration will then proceed with the new
preconditioner M . Moreover, since we compute the residual for each individual
column m k of M , it is easy to single out the most difficult columns
and concentrate on them to improve the convergence of the iterative solver.
This may prove useful in connection with flexible preconditioning, where the
preconditioner is adapted during the iterative process (see [17], [14], [3]).
Theoretical Properties of M
We shall now derive rigorous bounds on the spectrum of the preconditioned
matrix AM . Furthermore, we shall estimate the difference between M and
derive conditions that guarantee that M is nonsingular. Let M be
an approximate inverse of A obtained from the spai algorithm, and let m k
be its k-th column. We denote by r k the residual for every m k , and assume
that it satisfies
Theorem 3.1
1-k-n
fnumber of nonzero elements of r k
Then
We note that p is usually much smaller than n because A and M are
sparse. Thus, some of the bounds derived in this section are tighter than
those discussed in [4] and in (ch. 8, [1]). They directly apply to the computed
approximate inverse and the preconditioned system (2), because the
basic assumption (18) coincides with the stopping criterion used in the spai
algorithm.
Proof
Since
we immediately obtain the left inequality in (19). To derive the left inequality
in (20), we use the definition of the 2-norm
We get the right inequalities in (19) and (20) using (23) and
which holds both in the 2-norm and the Frobenius norm. As kr k k 1 is bounded
by immediately get (21).If we apply Gershgorin's theorem to AM , we see that all eigenvalues lie inside
a disk of radius centered at 1. Therefore, M is nonsingular if
We summarize this result as a corollary:
Corollary 3.1 If
As we generally do not know p in advance, (21) is only useful after having
computed M . Since p - n, we see that if p
must be nonsingular.
This gives a criterion for choosing " for a given n, although in practice it is
often too costly to run the algorithm with such a small ".
Even if A is symmetric, M will be nonsymmetric in general. It may then
be appropriate to use the symmetrized preconditioner
of M . To derive some estimates we begin with the following inequality, which
holds in the 2-norm and the Frobenius norm if A is symmetric:
Therefore, the symmetrized preconditioner (M +M T )=2 satisfies
In general, equation (26) is a very pessimistic estimate and not of much practical
use. This is because the spai algorithm does not take advantage of the
symmetry of A, and does not yield a symmetric approximate inverse. It is
easy, however, to reformulate the algorithm to compute only the lower triangular
part of M . This yields a symmetric preconditioner, but the algorithm
then loses its inherent parallelism. An interesting alternative would be to
compute a factorized sparse approximate inverse as in [11], but to leave the
sparsity open like in the spai algorithm.
The convergence of most iterative methods heavily depends on the distribution
of the eigenvalues or the singular values of the preconditioned matrix
[6]. Indeed, if most eigenvalues are clustered about 1 and only a few outliers
are present, the convergence will generally be very fast. Thus, it is crucial to
derive estimates on the spectrum of AM to determine the theoretical effectiveness
of the preconditioner. Such estimates are summarized in the following
two theorems.
Theorem 3.2 Let
1-k-n
fnumber of nonzero elements of r k
Then, the eigenvalues - k of AM are clustered at 1 and lie inside a circle of
radius
Proof
Let QRQ T be a Schur decomposition of AM \Gamma I. Then
the eigenvalues - k of AM are clustered at 1. Next, we use Gershgorin's
theorem and (21) to conclude that all - k lie inside a disk of radius p
centered at 1.
As an immediate consequence of (28), we see that minimizing kAM \Gamma Ik F also
reduces A's departure from normality kNk 2
is a Schur decomposition of AM .
We now derive bounds for the singular values and the condition number
of AM .
Theorem 3.3 The singular values of AM are clustered at 1 and lie inside
the interval [1 \Gamma
"). Furthermore, if
then the condition number of AM satisfies
s
Proof
Since the singular values of AM are the square roots of the eigenvalues of
We then bound (31) from above by
Next, we apply (28) to AMM T A T instead of AM , and conclude that the
singular values must be clustered at 1.Essential properties for the convergence of iterative methods are the clustering
of eigenvalues and singular values, the condition number, and the departure
from normality of the preconditioned linear system. In this section
we have shown that minimizing kAM \Gamma Ik in the Frobenius norm produces
a preconditioner M , which improves on all four points.
4 A first numerical example
We begin with a detailed study of orsirr2, the smallest among the orsx oil
reservoir simulation problems in the Harwell-Boeing matrix collection. Convergence
results for orsreg1 and orsirr1 are presented at the end of this section.
Here is a brief description of the problems:
bcg cgs bi-cgstab gmres(20) cgne

Table

1. Convergence results for orsirr2: unpreconditioned
preconditioned (0:2 - 0:6).
orsreg1: an oil reservoir simulation matrix for a 21 \Theta 21 \Theta 5 full grid, of size
entries.
orsirr1: an oil reservoir simulation matrix for an 21 \Theta 21 \Theta 5 irregular grid,
orsirr2: an oil reservoir simulation matrix for a 21 \Theta 21 \Theta 5 irregular grid,
Because of the rather small size of orsirr2, we can compute its true inverse
and all the eigenvalues and singular values of the preconditioned system. All
the numerical experiments presented in this section were computed in double
precision using a matlab implementation. The initial guess was always
and the stopping criterion
In table 1 we present the convergence results for a variety of iterative
methods and decreasing values of " . The algorithms for the various iterative
methods were obtained from [2]. The convergence results without preconditioning
are listed under I. We denote by cgne the cg algorithm applied
to the preconditioned normal equations AMM T A T y.
As we reduce " and kAM \Gamma Ik F , the number of iterations drops significantly
for all iterative methods. The results in table 1 exemplify the robustness of
the preconditioner with respect to ".

Table

2. Orsirr2: kAM \Gamma Ik and the condition number of AM for different
values of " .
In table 2 we compare different norms of AM \Gamma I as a function of ". As
we reduce ", both the 1-norm and the 2-norm remain nearly constant. Thus,
is a better solution of the minimization problem kAM \Gamma Ik 1 than
the computed M for all " - 0:2 . This is also true for the 2-norm as long
as " - 0:3, and indicates that neither the 1-norm nor the 2-norm is a good
measure of the proximity of AM to the identity, if the Frobenius norm is being
minimized. As we reduce " the singular values cluster at one, as demonstrated
in theorem 3.3. This is clearly shown in table 2, since the condition number
of AM is the ratio of the largest to the smallest singular value.
In figure 1 we verify the improvement in the clustering of the eigenvalues
predicted by theorem 3.2 as we reduce ". Next, we compare the approximate
bcg cgs bi-cgstab gmres(20) gmres(50)

Table

3. Convergence results for orsreg1 and orsirr1: unpreconditioned
preconditioned
-0.2
-0.2
Figure

1. Orsirr2: eigenvalues of AM with
0:2 (bottom).
inverse M with the true inverse A \Gamma1 . We compute A \Gamma1 and discard all entries
whose absolute value is less than or equal to 0.001 . Then, we compute the
approximate inverse M with 0:2 . It is quite striking how well the sparsity
patterns of both matrices agree with each other qualitatively in figure 2.
We conclude this section with the convergence results for both orsreg1
and orsirr1, given in table 3. The relative sparsity nz(M)=nz(A) was 0:94
for orsreg1 and 0:88 for orsirr1.
5 Numerical Experiments
In this section we consider a wide spectrum of problems coming from scientific
and industrial applications. We shall demonstrate the effectiveness
of the preconditioner for two standard but very different iterative methods:
bi-cgstab and gmres(m) with restart. We recall that the former requires
two matrix-vector multiplications per iteration, whereas the latter requires
only one matrix-vector multiply per iteration. In all numerical calculations

Figure

2. Orsirr2: sparsity pattern of entries in A \Gamma1 larger than 0.001 in
absolute value (left), and sparse approximate inverse M with
0:2 (right).
the initial guess was x unless specified, the stopping
criterion was as in (33). All the computations were done in double precision
fortran, partly on a Sparc10 sun station and partly on an Iris 4D/35 sgi
station. The fortran code provided in [2] was used for gmres.
We begin with the convergence results for the five shermanx black oil
simulators. This set consists of
sherman1: a black oil simulator, shale barrier,
of size
sherman2: a thermal simulation, steam injection, 6 \Theta 6 \Theta 5 grid, 5 unknowns,
sherman3: a black oil, impes simulation, 35 \Theta 11 \Theta 13 grid, 1 unknown,
sherman4: a black oil, impes simulation, 16 \Theta 23 \Theta 3 grid, 1 unknown,
sherman5: a fully implicit black oil simulator, 16 \Theta 23 \Theta 3 grid, 3 unknowns,
sh1 sh2 sh3 sh4 sh5

Table

4. Convergence results for shermanx: unpreconditioned (top), and
preconditioned (bottom).
The right-hand side was always provided. Table 4 shows that preconditioning
clearly improves the convergence for all considered problems. Here "max.
denotes the upper limit on the number of nonzero elements per
column in M . For sherman2, both bi-cgstab and gmres(20) reduced the
relative residual below 10 \Gamma5 after 4 and 7 steps respectively, but never reached
\Gamma8 . This may be due to the very large condition number 9:64 \Theta 10 11 of
sherman2, and could not be improved upon. A different implementation of
gmres [16] might further reduce the residual, but we did not pursue this
matter. In figure 3 we have displayed both the original matrix A and the
approximate inverse M for sherman2. This picture clearly shows why we
cannot simply set M equal to a banded matrix like in [10]. Indeed the sparsity
structures of A and M are totally different in this particular case, whereas
the number of nonzero entries in both matrices are comparable.
The next set of examples consists of the larger problems in the poresx and
collections:
pores2: a nonsymmetric matrix,
pores3: a nonsymmetric matrix,
saylor3: a nonsymmetric problem,
saylor4: a nonsymmetric problem,
Here the right-hand side was randomly chosen. Since pores2 did not generate
a very sparse approximate inverse and the iteration never reached the relative

Figure

3. Sherman2: original matrix A (left), and sparse approximate inverse
tolerance of 10 \Gamma8 , we opted for left instead of right preconditioning. Thus
in the case of pores2, we considered kMA \Gamma Ik F instead of (4). Yet we still
computed the exact residual of the original problem to check convergence.
This approach proved to be more efficient, because the rows of A \Gamma1 could be
approximated more effectively than its columns. In figure 4 we compare both
pores2 and its left approximate inverse. The number of nonzero entries is
similar, but the sparsity patterns are quite different. For pores2, gmres(20)
reduced the relative residual below 10 \Gamma6 after 106 iterations, but did not
improve any further. Again, a different implementation of gmres [16] might
mitigate this problem. In saylor3 we discovered in columns 988, 989, 998,
and 999 two independent 2 \Theta 2 singular submatrices. We ignored those four
columns in the computation of M , and simply replaced the two submatrices
by 2 \Theta 2 identity matrices. To guarantee the existence of a solution, we set the
right-hand side to A times a random vector. The orsx, shermanx, poresx, and
saylorx problems were all taken from the Harwell-Boeing matrix collection.
A comparative study of these problems for different iterative methods using
incomplete LU preconditioning can be found in [15].
In the final part of this section we shall consider several large problems.
The right-hand side was always provided, the initial guess x and the
tolerance as in (33). We start with four typical problems from Centric Engineering


Figure

4. Pores2: original matrix A (left), and sparse approximate inverse
pores 2 pores 3 saylor 3 saylor 4
bi-cgstab 78 118 104 285

Table

5. Convergence results for poresx and saylorx: unpreconditioned
(top), and preconditioned (middle & bottom). Left preconditioning
is used for pores2.

Table

6. Convergence results for Px: unpreconditioned (top), and preconditioned
(bottom).
P1: incompressible flow in a pressure driven pipe,
P2: incompressible flow in a pressure driven pipe,
P3: landing hydrofoil airplane FSE model,
P4: slosh tank model,
We remark that we have used the actual number of nonzero elements in
the Px matrices, since about 0.3% of the entries were zero in the original
data files. It is quite remarkable that although the Px matrices are rather
full, we obtain such a big improvement in convergence with much sparser
approximate inverses. Moreover, because of the periodic pattern in figure 5,
one could determine the sparsity pattern of the first few columns of M , and
then slide the pattern about the diagonal down to the lower right corner to
get the full sparsity structure.
To conclude this series of numerical experiments, we consider three problems
coming from an implicit 2-D Euler solver for an unstructured grid [19].
The matrices are of order
and correspond to the initial and later stages in the flow simulation. None of
the problems converged without preconditioning. The first problem T01 was
much easier, since the flow was still in its initial stage. It is also discussed
in [7]. Again, we see with T01 that we obtain a considerable improvement

Figure

5. P1: original matrix A (left), and sparse approximate inverse M
in convergence with a very sparse approximate inverse. We ran T01 with
0:4 and a maximum of 50 nonzero entries per column in M , which lead
to a relative sparsity of Both T25 and T50 ran with
0:3 and a maximum of 100 nonzero entries per column in M , which lead
to a relative sparsity of
Figure

6. Convergence history of Tx with bi-cgstab: relative residual vs.
number of iterations for T01(dash-dot), T25(dotted), T50(solid).
6 Conclusion
The spai algorithm computes a sparse approximate inverse M of a general
sparse matrix A. It is inherently parallel, since the columns of M are calculated
independently of one another. The matrix M gives valuable insight
into A \Gamma1 , and provides a measure on the proximity of M to A \Gamma1 . Instead of
imposing an a priori sparsity pattern upon M , we let the algorithm capture
automatically the relevant entries in the inverse. Thus, we minimize the number
of nonzero entries in M and concentrate the computational effort where
it is needed. The algorithm generates a robust and flexible preconditioner
for iterative solvers, as it gives full control over the sparsity and the quality
of M . We have shown both from a theoretical and a practical point of view
that the preconditioner generated by the spai algorithm is very effective in
improving the convergence of iterative solvers.
It is possible to minimize the total execution time for a particular problem
and architecture, by choosing an optimal M . It is clear that a very sparse
preconditioner is very cheap but may not lead to much improvement in con-
vergence, and that if M becomes too dense, it becomes too expensive to
compute. The optimal preconditioner lies somewhere between these two ex-
tremes, and is problem and architecture dependent. In a parallel environment
with very slow communication capabilities, such as a cluster of workstations,
it may be advantageous to compute a fairly dense approximate inverse. This
will increase the local floating-point intensive computations, and reduce the
number of iterations; hence, the amount of communications required at each
iteration by inner products and matrix-vector products will be diminished.
The implementation of this method on a parallel computer is straight-
forward, since all calculations are done independently. Yet, each processor
must have access to the data required to solve its subproblem. Although the
preconditioner is invariant under permutation, the ordering of the unknowns
can affect the amount of inter-processor communication involved in the computation
of M , and may be optimized for a particular application. This
approach will prove particularly effective when a linear system needs to be
solved repeatedly, such as in implicit time-marching schemes or the solution
of nonlinear equations.

Acknowledgements

We thank Horst Simon for providing us with the matrices and for his helpful
comments. This paper was written while the second author was visiting the
SCCM program at Stanford University. He would like to thank Gene Golub
for his kind hospitality.



--R

Cambridge University Press
Templates for the Solution of Linear Systems
Approximate inverse preconditioners for general sparse matrices
Approximate inverse preconditionings for sparse linear systems
Decay rates for inverses of band matrices
Iterative solution of linear systems
Implementation details of the coupled QMR algorithm
Numerical methods for solving linear least squares problems
Parallel preconditioning and approximate inverses on the Connection Machine
Parallel preconditioning and approximate inverses on the Connection Machine
Factorized sparse approximate inverse preconditionings
Factorized sparse approximate inverse (FSAI) preconditionings for solving 3D FE systems on massively parallel computers II
Sparse approximate inverse preconditionings for solving 3D CFD problems on massively parallel computers
A flexible inner-outer preconditioned GMRES algorithm
A comparative study of preconditioned Lanczos methods for nonsymmetric linear systems
Efficient high accuracy solutions with GM- RES(m)
GMRESR: a family of nested GMRES methods

Unstructured grid solvers on the iPSC/860
--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
Sangback Ma, Comparisons of the ILU(0), Point-SSOR, and SPAI Preconditioners on the CRAY-T3E for Nonsymmetric Sparse Linear Systems Arising from PDEs on Structured Grids, International Journal of High Performance Computing Applications, v.14 n.1, p.39-48, February  2000
Thomas Huckle, Factorized Sparse Approximate Inverses for Preconditioning, The Journal of Supercomputing, v.25 n.2, p.109-117, June
Ravindra Boojhawon , Muddun Bhuruth, Restarted simpler GMRES augmented with harmonic Ritz vectors, Future Generation Computer Systems, v.20 n.3, p.389-397, April 2004
J. Martnez , G. Larrazbal, Wavelet-based SPAI pre-conditioner using local dropping, Mathematics and Computers in Simulation, v.73 n.1, p.200-214, 6 November 2006
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
T. Tanaka , T. Nodera, Effectiveness of approximate inverse preconditioning by using the MR algorithm on an origin 2400, Proceedings of the third international conference on Engineering computational technology, p.115-116, September 04-06, 2002, Stirling, Scotland
Davod Khojasteh Salkuyeh , Faezeh Toutounian, BILUS: a block version of ILUS factorization, The Korean Journal of Computational & Applied Mathematics, v.15 n.1-2, p.299-312, May 2004
Oliver Brker , Oscar Chinellato , Roman Geus, Using Python for large scale linear algebra applications, Future Generation Computer Systems, v.21 n.6, p.969-979, June 2005
Y. Liang , J. Weston , M. Szularz, Generalized least-squares polynomial preconditioners for symmetric indefinite linear equations, Parallel Computing, v.28 n.2, p.323-341, February 2002
Olivier Besson, Band preconditioners: application to preconditioned conjugate gradient methods on parallel computers, Parallel numerical linear algebra, Nova Science Publishers, Inc., Commack, NY, 2001
George A. Gravvanis , Konstantinos M. Giannoutakis, Parallel preconditioned conjugate gradient square method based on normalized approximate inverses, Scientific Programming, v.13 n.2, p.79-91, April 2005
N. Guessous , O. Souhar, Multilevel block ILU preconditioner for sparse nonsymmetric M-matrices, Journal of Computational and Applied Mathematics, v.162 n.1, p.231-246, 1 January 2004
M. Bhuruth , M. K. Jain , A. Gopaul, Preconditioned iterative methods for the nine-point approximation to the convection---diffusion equation, Journal of Computational and Applied Mathematics, v.138 n.1, p.73-92, 1 January 2002
Kai Wang , Jun Zhang, Multigrid treatment and robustness enhancement for factored sparse approximate inverse preconditioning, Applied Numerical Mathematics, v.43 n.4, p.483-500, December 2002
Stephen T. Barnard , Luis M. Bernardo , Horst D. Simon, An MPI Implementation of the SPAI Preconditioner on the T3E, International Journal of High Performance Computing Applications, v.13 n.2, p.107-123, May       1999
Claus Koschinski, New methods for adapting and for approximating inverses as preconditioners, Applied Numerical Mathematics, v.41 n.1, p.179-218, April 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
Matthias Bollhfer , Volker Mehrmann, Some convergence estimates for algebraic multilevel preconditioners, Contemporary mathematics: theory and applications, American Mathematical Society, Boston, MA, 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
Prasanth B. Nair , Arindam Choudhury , Andy J. Keane, Some greedy learning algorithms for sparse regression and classification with mercer kernels, The Journal of Machine Learning Research, 3, 3/1/2003
George A. Gravvanis, On the Solution of Boundary Value Problems by Using Fast Generalized Approximate Inverse Banded Matrix Techniques, The Journal of Supercomputing, v.25 n.2, p.119-129, June
Michele Benzi , Miroslav Tma, A parallel solver for large-scale Markov chains, Applied Numerical Mathematics, v.41 n.1, p.135-153, April 2002
Keita Teranishi , Padma Raghavan , Esmond Ng, A new data-mapping scheme for latency-tolerant distributed sparse triangular solution, Proceedings of the 2002 ACM/IEEE conference on Supercomputing, p.1-11, November 16, 2002, Baltimore, Maryland
G. A. Gravvanis , K. M. Giannoutakis , M. P. Bekakos , O. B. Efremides, Parallel and systolic solution of normalized explicit approximate inverse preconditioning, The Journal of Supercomputing, v.30 n.2, p.77-96, November 2004
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
F. F. Hernndez , J. E. Castillo , G. A. Larrazbal, Large sparse linear systems arising from mimetic discretization, Computers & Mathematics with Applications, v.53 n.1, p.1-11, January, 2007
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
J. M. Dennis , H. M. Tufo, Scaling climate simulation applications on the IBM Blue Gene/L system, IBM Journal of Research and Development, v.52 n.1, p.117-126, January 2008
Michele Benzi, Preconditioning techniques for large linear systems: a survey, Journal of Computational Physics, v.182 n.2, p.418-477, November 2002
parallel algorithms;approximate inverses;sparse linear systems;sparse matrices;iterative methods;preconditioning