Sparse approximate inverse smoothers for geometric and algebraic multigrid.
Sparse approximate inverses are considered as smoothers for geometric and algebraic multigrid methods. They are based on the SPAI-Algorithm [MJ. Grote, T. Huckle, SIAM J. Sci. Comput. which constructs a sparse approximate inverse M of a matrix A, by minimizing I - MA in the Frobenius norm. This leads to a new hierarchy of inherently parallel smoothers: SPAI-0, SPAI-1, and SPAI(). For geometric multigrid, the performance of SPAI-1 is usually comparable to that of Gauss-Seidel smoothing. In more difficult situations, where neither Gauss-Seidel nor the simpler SPAI-0 or SPAI-1 smoothers are adequate, further reduction of  automatically improves the SPAI() smoother where needed. When combined with an algebraic coarsening strategy [J.W. Ruge, K. Stben, in: S.F. McCormick (Ed.), Multigrid Methods, SIAM, 1987, pp. 73-130] the resulting method yields a robust, parallel, and algebraic multigrid iteration, easily adjusted even by the non-expert. Numerical examples demonstrate the usefulness of SPAI smoothers, both in a sequential and a parallel environment.Essential advantages of the SPAI-smoothers are: improved robustness, inherent parallelism, ordering independence, and possible local adaptivity.
Introduction
Multigrid methods rely on the subtle interplay of smoothing and coarse grid
correction. Only their careful combination yields an e-cient multigrid solver
for large linear systems, resulting from the discretization of partial dierential
equations [7,14,15,27]. Standard smoothers for multigrid usually consist of a
few steps of a basic iterative method. Here we shall consider smoothers that are
based on sparse approximate inverses. Hence, starting from the linear system
we let M denote a sparse approximation of A 1 . The corresponding basic
iterative method is
Since the approximate inverse M is known explicitly, each iteration step requires
only one additional M  v matrix-vector multiply; therefore, it is easy
to parallelize and cheap to evaluate, because M is sparse.
Recently, various algorithms have been proposed, all of which attempt to compute
directly a sparse approximate inverse of A [5,9,17]. For a comparative
study of various approximate inverse preconditioners we refer to Benzi and
Tuma [6]. Approximate inverse techniques are also gaining in importance as
smoothers for multigrid methods. First introduced by Benson and Frederickson
[3,4], they were shown to be eective on various di-cult elliptic problems
on unstructured grids by Tang and Wan [25]. Advantages of sparse approximate
inverse smoothers over classical smoothers, such as damped Jacobi,
Gauss-Seidel, or ILU, are inherent parallelism, possible local adaptivity, and
improved robustness.
Here we shall consider sparse approximate inverse (SPAI) smoothers based on
the SPAI-Algorithm by Grote and Huckle [13]. The SPAI-Algorithm computes
an approximate inverse M explicitly by minimizing I MA in the Frobenius
norm. Both the computation of M and its application as a smoother are inherently
parallel. Since an eective sparsity pattern of M is in general unknown a
priori, the SPAI-Algorithm attempts to determine the most promising entries
dynamically. This strategy has proved eective in generating preconditioners
for many di-cult and ill-conditioned problems [1,13,24]. Moreover, it provides
the means for adjusting the smoother locally and automatically, if necessary.
Nevertheless, by choosing an a priori sparsity pattern for M , the computational
cost can be greatly reduced. Possible choices include powers of A or
A > A, as suggested by Huckle [16] or Chow [10]. Hence we shall investigate the
following hierarchy of sparse approximate inverse smoothers: SPAI-0, SPAI-1,
and SPAI("). For SPAI-0 and SPAI-1 the sparsity pattern of M is xed: M is
diagonal for SPAI-0, whereas for SPAI-1 the sparsity pattern of M is that of
A. For SPAI(") the sparsity pattern of M is determined automatically by the
SPAI-Algorithm ([13]); the parameter " controls the accuracy and the amount
of ll-in of M .
As structured geometric grids are di-cult to use with complex geometries,
application code designers often turn to very large unstructured grids. Yet the
lack of a natural grid hierarchy prevents the use of standard geometric multi-
grid. In this context, algebraic multigrid (AMG) is often seen as the most
promising method for solving large-scale problems. The original AMG algo-
rithm, rst introduced in the 1980's by Ruge and Stuben [22], uses the (simple)
Gauss-Seidel iteration as a smoother, but determines the coarse "grid" space
in a sophisticated way to improve robustness of the method. But if the iteration
fails to converge, there is no automatic way to improve on the smoother.
As an alternative, we investigate the usefulness of smoothers based on SParse
Approximate Inverses (SPAI). Not only inherently parallel, their performance
can also easily be adjusted, even by a non-expert. Thus we aim for a more
general and inherently parallel algebraic multigrid method.
In Section 2 we brie
y review the SPAI-Algorithm and show how sparse
approximate inverses are used as smoothers within a multigrid iteration. A
heuristic Green's function interpretation underpins their eectiveness as smoo-
thers. Rigorous results on the smoothing property of approximate inverses were
proved in [8]; they are summarized in Section 2.4. Next, we present in Section
3 a detailed description of the algebraic coarsening strategy used ([22]),
together with key components of the algorithm for e-cient implementation.
Finally, in Section 4, we compare the performance of SPAI smoothing to that
of Gauss-Seidel smoothing on various test problems, either within a geometric
or an algebraic multigrid setting.
smoothing
2.1 Classical smoothers
Consider a sequence of nested grids, T On the nest mesh,
we wish to solve the n  n linear system
by a multigrid method { for further details on multigrid see Hackbusch [14]
and ([15], Section 10) or Wesseling [27]. A multigrid iteration results from the
recursive application of a two-grid method. A two-grid method consists of  1
pre-smoothing steps on level ', a coarse grid correction on level ' 1, and  2
post-smoothing steps again on level '. This leads to the iteration
for the error e (m)
and r are prolongation and restriction
operators, respectively, between T ' and T ' 1 , while S ' denotes the smoother. If
nested nite element spaces with Galerkin discretization are used, the Galerkin
product representation holds:
Otherwise, one can still use (4) to dene the coarse-grid problem for given r
and p.
We shall always use x (0)
our initial guess. The multigrid (V-cycle) iteration
proceeds until the relative residual drops below a prescribed tolerance,
kb A ' x (m)
kbk < tol: (5)
Then we calculate the average rate of convergence
1=m
The expected multigrid convergence behavior is achieved if the number of
multigrid iterations, m, necessary to achieve a xed tolerance, is essentially
independent of the number of grid levels '.
Typically, the smoother has the form
where W ' approximates A ' and is cheap to invert | of course, W 1
' is never
computed explicitly. Let with D the diagonal, L the lower
triangular part, and U the upper triangular part of A. Then, damped Jacobi
smoothing corresponds to
whereas Gauss-Seidel smoothing corresponds to
In (8) the parameter ! is chosen to maximize the reduction of the high frequency
components of the error. The optimal value, !  , is problem dependent
and usually unknown a priori. Although Gauss-Seidel typically leads to faster
convergence, it is more di-cult to implement in parallel, because each smoothing
step in (9) requires the solution of a lower triangular system. If neither
Jacobi nor Gauss-Seidel smoothing lead to satisfactory convergence, one can
either resort to more sophisticated (matrix dependent) prolongation and restriction
operators ([29]) or to more robust smoothers, based on incomplete LU
(ILU) factorizations of A ' ([28]). Unfortunately, ILU-smoothing is inherently
sequential and therefore di-cult to implement in parallel. It is also di-cult
to improve locally, say near the boundary or a singularity, without aecting
the ll-in everywhere in the LU factors.
2.2 SPAI-smoothers
As an alternative to inverting W ' in (7), we propose to compute explicitly a
sparse approximate inverse M of A, and to use it for smoothing { we drop the
grid index ' to simplify the notation. This yields the SPAI-smoother,
where M is computed by minimizing kI MAk in the Frobenius norm for
a given sparsity pattern. In contrast to W 1
' in (7), the matrix M in (10) is
computed explicitly. Therefore, the application of S SPAI requires only matrix-vector
multiplications, which are easy to parallelize; it does not require the
solution of any upper or lower triangular systems. Moreover, the Frobenius
norm naturally leads to inherent parallelism because the rows m >
of M can
be computed independently of one another. Indeed, since
the solution of (11) separates into the n independent least-squares problems
for the m >
Here e k denotes the k-th unit vector. Because A and M are sparse, these
least-squares problems are small.
Since an eective sparsity pattern of M is unknown a priori, the original SPAI-
Algorithm in [13] begins with a diagonal pattern. It then augments progressively
the sparsity pattern of M to further reduce each residual r
Each additional reduction of the 2-norm of r k involves two steps. First, the
algorithm identies a set of potential new candidates, based on the sparsity of
A and the current (sparse) residual r k . Second, the algorithm selects the most
protable entries, usually less than ve entries, by computing for each candidate
a cheap upper bound on the reduction in kr k k 2 . Once the new entries
have been selected and added to m k , the (small) least-squares problem (12) is
solved again with the augmented set of indices. The algorithm proceeds until
each row m >
of M satises
where " is a tolerance set by the user; it controls the ll-in and the quality of
the preconditioner M . A larger value of " leads to a sparser and less expensive
approximate inverse, but also to a less eective smoother with a higher number
of multigrid cycles. A lower value of " usually reduces the number of cycles,
but the cost of computing may become prohibitive; moreover,
a denser M results in a higher cost per smoothing step. The optimal value of
" minimizes the total time; it depends on the problem, the discretization, the
desired accuracy, and the computer architecture. Further details about the
original SPAI-Algorithm can be found in [13].
In addition to SPAI("), we shall also consider the following two greatly sim-
plied SPAI-smoothers with xed sparsity patterns: SPAI-0, where M is di-
agonal, and SPAI-1, where the sparsity pattern of M is that of A. Both solve
the least-squares problem (12), and thus minimize kI MAk in the Frobenius
norm for the sparsity pattern chosen a priori. This eliminates the search
for an eective sparsity pattern of M , and thus greatly reduces the cost of
computing the approximate inverse. The SPAI-1 smoother coincides with the
smoother of Tang and Wan [25].
For diagonal and can be calculated directly. It is
simply given by
a kk
where a >
k is the k-th row of A { note that M is always well-dened if A is
nonsingular. In contrast to damped Jacobi, SPAI-0 is parameter-free.
To summarize, we shall consider the following hierarchy of SPAI-smoothers,
which all minimize kI MAk in the Frobenius norm for a certain sparsity
pattern of M .
diagonal and given by (14).
SPAI-1: The sparsity pattern of M is that of A.
SPAI("): The sparsity pattern of M is determined automatically via the SPAI-
Algorithm [13]. Each row m >
of M satises (13) for a given ".
Any of these approximate inverses leads to the smoothing step
We have found that in many situations, SPAI-0 and SPAI-1 yield ample
smoothing. However, the added
exibility in providing an automatic criterion
for improving the smoother via the SPAI-Algorithm remains very useful.
Indeed, both SPAI-0 and SPAI-1 can be used as initial guess for SPAI("), and
thus be locally improved upon where needed by reducing ". For matrices with
inherent (small) block structure, typical from the discretization of systems of
partial dierential equations, the Block-SPAI-Algorithm ([2]) greatly reduces
the cost of computing M .
2.3 Green's function interpretation
Why do approximate inverses yield eective smoothers for problems which
come from partial dierential equations? As the mesh parameter h tends to
zero, the solution of the linear system,
A h u
tends to the solution of the underlying dierential equation,
with appropriate boundary conditions. Here the matrix A h corresponds to a
discrete version of the dierential operator L. Let y h
k denote the k-th row of
solves the linear system
with e h
k the k-th unit vector. As h ! 0, y h
k tends to the Green's function
Here L  denotes the adjoint dierential operator and -(x x k ) the \delta-
centered about x k . To exhibit the correspondence between y h
k and
we recall that L  is formally dened by the identity
for all u; v in appropriate function spaces. Equation (20) is the continuous
counterpart to the relation
From (17), (19) and (20) we conclude that
Similarly, the combination of (16), (18), and (21) leads to the discrete counterpart
of (22),
Comparison of (22) and (23) shows that y h
k corresponds to to G(x; x k ) as
The k-th row, (m h
of the approximate inverse, M h , solves (12), or equiva-
lently
min
Hence m h
k approximates the k-th column of A > , that is y h
k in (18), in the
(discrete) 2-norm for a xed sparsity pattern of m h
k . The nonzero entries of
usually lie in a neighborhood of x k : they correspond to mesh points x j close
to x k . Therefore, after an appropriate scaling in inverse powers of h, we see
that m h
approximates G(x; x k ) locally in the (continuous) L
1268 3 For
(partial) dierential operators, G(x; x k ) typically is singular at x k and decays
smoothly, but not necessarily rapidly, with increasing distance jx x k j. Clearly
the slower the decay, the denser M h must be to approximate well A 1
deciency of sparse approximate inverse preconditioners was also pointed out
by Tang [24]. At the same time, however, it suggests that sparse approximate
inverses, obtained by the minimization of kI MAk in the Frobenius norm,
naturally yield smoothers for multigrid. Indeed to be eective, a preconditioner
must approximate uniformly over the entire spectrum of L. In contrast,
an eective smoother only needs to capture the high-frequency behavior of
. Yet this high-frequency behavior corresponds to the singular, local
behavior of G(x; x k ), precisely that which is approximated by m h
k .
To illustrate this fact, we consider the standard ve-point stencil of the discrete
Laplacian on a 15  15 grid. In Figure 1 on the following page we compare
A 1 with the Gauss-Seidel approximate inverse, (L and two explicit
approximate inverses, SPAI-1 and SPAI(0:2). We recall that Gauss-Seidel, a
poor preconditioner for this problem, remains an excellent smoother, because
it captures the high-frequency behavior of A 1 . Similarly, SPAI-1 and SPAI(")
yield local operators with, as we shall see, good smoothing property. Despite
the resemblance between the Gauss-Seidel and the SPAI approximate inverses,
we note the one-sidedness of the former, in contrast to the symmetry of the
latter.
Gauss-Seidel,
SPAI (0.2)
Fig. 1. Row 112 of the following operators: A 1 (top left), the Gauss-Seidel inverse
computed with SPAI-1 (bottom left), and M computed
with SPAI(0:2) (bottom right).
2.4 Theoretical Properties
In contrast to the heuristic interpretation of the previous section, we shall
now summarize some rigorous results ([8]) on the smoothing property of the
simplest smoother: SPAI-0.
Multigrid convergence theory rests on two fundamental conditions: the smoothing
property ([15], Denition 10.6.3):
any function with lim
and the approximation property ([15], Section 10.6.3). In general, the smoothing
and approximation properties together imply convergence of the two-grid
method and of the multigrid W-cycle, with a contraction number independent
of the level number '. Moreover, for symmetric positive denite prob-
lems, both conditions also imply multigrid V-cycle convergence independent
of ' { see Hackbusch ([15], Sect. 10.6) for details. The approximation property
is independent of the smoother, S ' ; it depends only on the discretization
the prolongation operator p, and the restriction operator r. In
[15] the approximation property is shown to hold for a large class of discrete
elliptic boundary value problems. For symmetric positive denite problems
the smoothing property usually holds for classical smoothers, such as damped
Jacobi, (symmetric) Gauss-Seidel, and incomplete Cholesky.
In [8] the smoothing property (25) was shown to hold for the SPAI-0 smoother
under reasonable assumptions on the matrix A. More precisely, for A symmetric
and positive denite, the SPAI-0 smoother satises the smoothing property,
either if A is weakly diagonally dominant, or if A has at most seven nonzero
o-diagonal entries per row.
Furthermore, the two diagonal smoothers SPAI-0 and damped Jacobi, with
optimal relaxation parameter !  , lead to identical smoothers for the discrete
Laplacian with periodic boundary conditions in any space dimension
[8]. In this special situation, the parameter-free SPAI-0 smoother automatically
yields a scaling of diag(A), which minimizes the smoothing factor; in
that sense it is optimal. In more general situations, however, both smoothers
dier because of boundary conditions, even with constant coe-cients on an
equispaced mesh. Comparison of these two diagonal smoothers via numerical
experiments showed that SPAI-0 is an attractive alternative to damped Jacobi
[8]. Indeed, SPAI-0 is parameter-free and typically leads to slightly better
convergence rates than damped Jacobi.
3 Algebraic Multigrid
Multigrid (MG) methods are sensitive to the subtle interplay between smoothing
and coarse-grid correction. When a standard geometric multigrid method
is applied to di-cult problems, say with strong anisotropy, this interplay is
disturbed because the error is no longer smoothed equally well in all direc-
tions. Although manual intervention and selection of coarse grids can sometimes
overcome this di-culty, it remains cumbersome to apply in practice to
unstructured grids and complex geometry. In contrast, an algebraic multigrid
approach compensates for the decient smoothing by a sophisticated
choice of the coarser grids and the interpolation operators, which is only based
on the matrix A ' . Many AMG variants exist, which dier in the coarsening
strategy or the interpolation used { an introduction to various AMG methods
can be found in ([26]).
Following Ruge and Stuben [22], we now describe the algebraic coarsening
strategy and interpolation operators, which we shall combine with the SPAI
smoothers from Section 2.2 and use for our numerical experiments.
3.1 Coarsening strategy1020301020300.20.6
x
y
Anisotropic stencil:h 26 6 6 6 6 4
Fig. 2. The error after ve Gauss-Seidel smoothing steps for the problem described
in Section 4.2 on a with 0:01. The smooth error component is aligned with the
anisotropy, which can be read from the stencils.
The fundamental principle underlying the coarsening strategy is based on the
observation that interpolation should only be performed along smooth error
components. For symmetric M-matrices, the error is smoothed well along large
(negative) o-diagonal entries in the matrix A ([23]). Therefore, at each grid
point p, we may identify among neighboring points q good candidates for
interpolation, by comparing the magnitude of the corresponding entries a pq .
This leads to the following relations between the point p and its neighbors q
in the connectivity graph of the matrix A:
Condition Notation Interpretation
a pq   max apr <0 ja pr j p ( q p (strongly) depends on q
and a pq 6= 0 or
q (strongly) in
uences p
a pq <  weakly depends on q
and a pq 6= 0 or
q weakly in
uences p
The parameter  controls the threshold, which discriminates between strong
and weak connections; typically 0:25. With this denition, all positive
o-diagonal entries are necessarily weak. The relations p ( q and p q are
symmetric only if A is symmetric.
Next, we dene the set of dependencies of a point p as
the set of in
uences of a point p as
I
and the set of weak dependencies of a point p as
On every level, the coarsening strategy must divide P , the set of all points
on that level, into two disjoint sets: C, the \coarse points", also present on
the coarser level, and F , the \ne points", which are absent on the coarser
level. The choice of C and F induces the C/F{splitting of
Coarse grid correction heavily depends on accurate interpolation. Accurate
interpolation is guaranteed if every F point is surrounded by su-ciently many
strongly dependent C points. A typical conguration is shown in Figure 3.
q3
strong
dependency
dependency
C point
F point
Strong dependencies are indicated with
solid arrows, while weak dependencies are
represented by dashed arrows. C points
are represented by solid circles, whereas F
points are represented by dashed circles.
Hence all the strong dependencies of point
are C points; therefore q 2 and
q 4 are good candidates for interpolating p.
Fig. 3. Ideal coarsening conguration for interpolation
The coarsening algorithm attempts to determine a C/F{splitting, which maximizes
the F {to{C dependency for all F points (Coarsening Goal 1), with
a minimal set C (Coarsening Goal 2). It is important to strike a good balance
between these two con
icting goals, as the overall computational eort
depends not only on the convergence rate, but also on the amount of work
per multigrid cycle. Clearly the optimal C=F -splitting minimizes total execution
time. However, since the convergence rate is generally unpredictable, the
coarsening algorithm merely attempts to meet Coarsening Goals 1 and 2 in
a heuristic fashion. In doing so, its complexity must not exceed O(n log n) to
retain the overall complexity of the multigrid iteration.
3.2 Coarse grid selection: a greedy heuristic
To split P into C and F , every step of a greedy heuristic moves the most
promising candidate from P into C, while forcing neighboring points into F .
This procedure is repeated until all points are distributed. If every step requires
at most O(log n) operations, and the complexity of all other computations does
not exceed O(n log n), the desired overall complexity of O(n log n) is reached.
The greedy heuristic described in [23] is based on the following two principles,
which correspond to Coarsening Goals 1 and 2:
(1) The most promising candidate, p, for becoming a C point, is that with
the highest number of in
uences jI p j. Then all in
uences of p are added
to F . This choice supports Coarsening Goal 1 because all F points will
eventually have at least one strong C dependency.
(2) To keep the number of C points low (Coarsening Goal 2), the algorithm
should prefer C points near recently chosen F points; hence, these in
u-
ences are given a higher priority.
Starting with all points as \undecided points", that is the algorithm
proceeds by selecting from U the most promising C point with highest priority.
The priority of any point p is dened by
Equation re
ects the preference in choosing the next C point for a point
which in
uences many previously selected F points. The key advantage
of (26) is the possibility to update the priority locally and in O(1) time, which
results in the desired overall complexity of O(n log n). We now summarize the
Coarse Grid Selection algorithm:
Algorithm 1 Coarse Grid Selection
procedure
for all
set Priority(p) := jI p j
end for all
U :=
while U
(1) select p 2 U with maximal P riority(p)
for all q 2 D p (all dependencies of p)
end for all
for all q 2 I p (all in
uences of p)
for all r 2 D q (all dependencies of q)
end for all
end for all
while
procedure
To implement steps (1), (2), and (3) e-ciently in O(1) time, we maintain a
list Q of all points sorted by priority, together with the list I of point indices
of Q. Moreover, a list of boundaries B of all priorities occurring in Q enables
the immediate update of the sorted list Q. Figure 4 shows a possible segment
of the lists Q, B, and I.
5 611 15index of Q
position in Q
priority(p)
position in Q3141312
I
index of B
Fig. 4. The three lists Q, I, and B enable the e-cient implementation of the coarsening
algorithm.
During the set{up phase of the coarsening algorithm, B is computed and
sorted by priority. Step (1) simply chooses the last element of Q. Steps (2)
and (3) are implemented by exchanging a point, whose priority must be either
incremented or decremented, with its left{ or rightmost neighbor in Q with
that same priority; then its priority is adjusted, while Q remains sorted. Both
B and I are updated accordingly. Following a suggestion of K. Stuben, we
shall skip the second pass of the original coarsening algorithm in [22], which
enforces even stronger F {to{C dependency, because of the high computational
cost involved.
3.3 Interpolation
The grid function u h , dened on the ner grid, is interpolated from the grid
function uH , dened on the coarser grid C, as follows:
Hence, values at C points are simply transfered from the coarser level, whereas
values at F points are interpolated from C neighbors. The four dierent dependencies
possible between any F point and its neighbors are shown in Figure 5.
For \standard interpolation" (see [23]), the choice of the weights, w pq , for
interpolating p, is based on the equation
a pp e
a pq e
Indeed, if Ae ' 0, the smoothing eect is minimal, and the error e is declared
\algebraically smooth" { see [23] for details. Clearly, we cannot interpolate p
from surrounding F points, whereas weakly dependent C points are not included
either, because of the rough nature of the error in that direction. Thus
connections (q 1 and q 2 in Figure 5) are always ignored in the interpolation
and the corresponding interpolation weights set to zero (w p;q 1
cancellation of the weak dependencies, the neglected entries of the weakly
dependent neighbors are added to the diagonal. Hence equation (28) becomes
~ a pp e
a pq e
a
Strong C dependencies, such as q 3 in Figure 5, cause no di-culty because the
value of uH is available at that coarse grid location. Division of (29) by ~ a pp
yields the weight
a pq
~ a pp
However, strong F dependencies, such as q 4 in Figure 5, are not available for
interpolation and must rst be interpolated from C points, on which they
strongly depend. To do so, we replace a qq by ~
a qq for every q 2 D q \ F , with
~ a
For every point r 2 D q this yields the weight
a pq
~ a pp
a qr
~ a qq
If a a point q is both a direct and an indirect neighbor of p, so that both (30)
and (32) apply, the two weights are calculated separately and then added to
each other.
q4
strong
dependency
dependency
C point
F point
q3
indirect interpolation
direct interpolation
ignored ignored
Fig. 5. The four dierent dependencies possible between p and its neighbors.
The algorithm described above determines the coarse grid levels only on the
basis of A, and not on that of the approximate inverse M . In fact, the information
contained in M can be used to determine coarse grid levels and
interpolation weights, as suggested by Meurant [19,20].
3.4 Measuring computational costs and memory requirements
When comparing the performance of various smoothers, we cannot limit ourselves
to comparing the number of multigrid iterations, but also need to estimate
the additional amount of work due to the smoother. To do so, we
calculate the total density ratio,  M , of nonzero entries in M to those in A on
all grid levels, 1  i  ', where smoothing is applied:
The additional amount of work due to the smoother is proportional to  M .
While rapidly reducing the number of points from one level to the next, the
matrices A i must also remain reasonably sparse, as measured by
For instance, as Galerkin coarse grid approximation enlarges the standard ve{
point stencil on the nest grid to nine{point stencils on subsequent levels, the
resulting value of  A for geometric multigrid is about 1.6. If semi{coarsening
together with one-dimensional interpolation is used,  A increases up to two.
All the results presented in the following section were computed with a MATLAB
implementation. We shall evaluate the e-ciency of the various approaches
by comparing their respective values for  M and  A .
4 Numerical results
To illustrate the usefulness and versatility of SPAI smoothing, we shall now
consider various standard test problems. In all cases, the dierential equation
considered is discretized on the nest level with standard nite dierences
on an equispaced mesh. For geometric multigrid, we use a regularly rened
sequence of equispaced grids, with a single unknown remaining at the center
of the domain. For algebraic multigrid, the coarser levels are obtained by the
Coarse Grid Selection Algorithm described in Section 3.2. With
the denition of strong dependency from section 3.1, the algorithm proceeds
until the number of grid points drops below twenty. The coarse grid operators
are obtained via the Galerkin product formula (4), with . For geometric
multigrid, p correspond to standard linear interpolation, whereas for AMG p
is obtained as described in Section 3.3. We use a multigrid V-cycle iteration,
with two pre- and two post-smoothing steps 2). The multigrid
iteration proceeds until the relative residual satises the prescribed tolerance
in (5), with
4.1 Rotating
ow problem
We rst consider the convection{diusion problem,
in (0; 1)(0; 1), with u(x; on the boundary. Here u represents any scalar
quantity advected by the rotating
ow eld. For convection dominated
<< h, the linear systems cease to be symmetric and positive denite, so
that these problems lie outside of classical multigrid theory. We use centered
second-order nite dierences for the diusion, but discretize the convection
with rst-order upwinding to ensure numerical stability.

Table
Geometric MG convergence rates for the rotating
ow problem on a 128128 grid,
for dierent values of . The symbol y indicates that the multigrid iteration diverges.
Smoother
Gauss-Seidel SPAI-0 SPAI-1 SPAI(0.3) SPAI(0.2)

Table

convergence rates q obtained with standard MG. All
smoothers yield acceptable convergence rates in the diusion dominated case,
with . For however, the multigrid iteration diverges with
Gauss{Seidel or SPAI-0 smoothing. In contrast, the SPAI-1 smoother still
yields a convergent method. The use of SPAI(0.3) smoothing accelerates convergence
even further, while  M increases up to 1.4 only.
As we reduce the diusion even further down to only the SPAI(0.2)
smoother yields a convergent iteration. Although the resulting value of  M
is quite high, the construction of SPAI(0.2) remains parallel and fully auto-
matic. We remark that symmetric Gauss-Seidel smoothing ([27]) leads to a
convergent multigrid iteration, yet this approach does not generalize easily to
unstructured grids.
Parallel results
Since the SPAI-1 smoother is inherently parallel, it is straightforward to apply
within a parallel version of geometric MG. The data is distributed among
processors via domain decomposition, which is well{known to work e-ciently
for a number of multigrid applications ([18]). The platform we shall use is
the ETH{Beowulf cluster, which consists of 192 dual CPU Pentium III (500
processors. All nodes are connected via a
100 MBit/s and 1 GB/s switched network, while communication is done with
MPI.
We now apply our parallel multigrid implementation to the rotating Flow
Problem (35) with On 128 nodes, the total execution{time is 156
seconds on the 40964096 grid. The time includes the set{up for the construction
of the SPAI-1 smoother, which requires the solution of about sixteen
million small (259) and independent least{squares problems. As shown
in table 1, the use of a coarsest level, which consists only of a single mesh
point, leads to a divergent multigrid iteration for increasing
the resolution of the coarsest level up to 3232 mesh points, one obtains a
convergent multigrid iteration.

Table
Scalability of parallel MG using SPAI-1. The problem size and the number of processors
is increased by a factor of 4, while total time increases by 30% only.
Gridsize 512512  4 10231023
Number of processors
Total time (sec) 20 26
To obtain good speed{up with a parallel MG code, it is important to perform
coarse grid agglomeration (see [21]) because of the loss of e-ciency on coarser
grid levels. Although we have not implemented such an agglomeration strategy,
our computations scale reasonably well as long as the problem size matches
the size of the parallel architecture { see Table 2.
Rotating
ow problem: algebraic
coarsening is clearly aligned
with the
ow direction. Larger dots
correspond to C points on coarser
levels.
(b) Locally anisotropic diusion:
semi-coarsening is apparent in the
center of the domain.
Fig. 6. Examples of algebraic coarsening for the two model problems considered.
AMG results
None of these approaches, however, is entirely satisfactory for vanishing vis-
cosity. To overcome the lack of robustness for small , we now apply the
algebraic coarsening strategy described in Section 3. Figure 6(a) displays the
coarse levels selected by the algorithm. In Table 3 both SPAI-0 and SPAI-1
yield convergence without any particular tuning. 0:5, the SPAI(")

Table
AMG convergence results for the rotating
ow problem for varying  on a 128128
grid.
A 2:8 3:4 4:2
Smoother q M q M q M
Gauss{Seidel 0.14 | 0.38 | 0.81 |

Table
AMG convergence rates for the rotating
ow problem with
Smoother
Gauss{Seidel SPAI-0 SPAI-1 SPAI(0.5)
Gridsize A q M q M q M q M
128128 4.2 0.81 | 0.36 (0.1) 0.21 (1.0) 0.27 (0.4)
4.3 0.96 | 0.38 (0.1) 0.24 (1.0) 0.34 (0.4)
smoother yields a compromise between the SPAI-0 and SPAI-1 smoothers:
both the storage requirement and the convergence rate lie between those obtained
with the xed sparsity patterns of SPAI-0 and SPAI-1. Lower values
of " reduce the convergence rate even further. The poor convergence rates
obtained with Gauss-Seidel could probably be improved, either by smoothing
C points before F points ([23]) or by using symmetric Gauss{Seidel.
The results in Table 4 demonstrate the robustness of SPAI smoothing. In-
deed, as  ! 0 all convergence rates obtained with the combined SPAI-AMG
approach remain bounded.
4.2 Locally anisotropic diusion
In this section we consider the locally anisotropic problem,

with u(x; on the boundary. The diusion coe-cient (x;
except inside the square [1=4; 3=4]  [1=4; 3=4], where (x; y)   is
constant. In Table 5 we observe that geometric multigrid has di-culties for
small values . Because of the unidirectional smoothing of the error, aligned
with the strong anisotropy, standard (isotropic) interpolation fails.

Table
Locally anisotropic diusion: geometric MG convergence rates q for varying  on a
128  128 grid.
Smoother q M q M q M
Gauss{Seidel
AMG results
AMG overcomes these di-culties by performing automatic semi{coarsening
and operator dependent interpolation only in the direction of strong couplings,
which correspond to smooth error components. It is well{known (e.g. [23])
that AMG solves such problems with little di-culty. The results in Table 6
verify this fact for acceptable densities  A and  M . Both
densities could be lowered even further by dropping the smallest entries in the
interpolation operators ([23]); we do not consider such truncated grid transfer
operators here.

Table
AMG convergence results for locally anisotropic diusion on a 128  128 grid. Note
that q, A , and M remain bounded as  ! 0.
A 2.84 2.94 2.94
Smoother q M q M q M
Gauss{Seidel 0.14 | 0.18 | 0.18 |
The convergence rates obtained with Gauss-Seidel and SPAI-1 are comparable
and both below 0.2, while SPAI-0 results in slightly slower convergence.
Overall the SPAI-1 smoother is the most e-cient smoother for this problem.
Although further reduction of " results in even faster convergence, the approximate
inverses become too dense and thus too expensive. Again, the results in

Tables

6 and 7 demonstrate robust multigrid behavior, either as h ! 0, or as

Table
AMG convergence rates q for locally anisotropic diusion, with . Note that
q; A , and M remain bounded as
Gridsize 6464 128128 256256
A 2.89 2.94 2.95
Gauss-Seidel 0.12 | 0.18 | 0.22 |
Concluding remarks
Our results show that sparse approximate inverses, based on the minimization
of the Frobenius norm, provide an attractive alternative to classical Jacobi
or Gauss-Seidel smoothing. The simpler smoothers, SPAI-0 and SPAI-1, often
provide ample smoothing, comparable to damped Jacobi or Gauss-Seidel. Nev-
ertheless, situations such as convection dominated rotating
ow, where SPAI-1
leads to a convergent multigrid iteration, unlike Gauss-Seidel, demonstrate the
improved robustness. Our implementation of geometric multigrid combined
with SPAI-1 smoothing enables the fast solution of very large convection-
diusion problems on massively parallel architectures. By incorporating the
SPAI smoothers into AMG ([22]), we obtain a
exible, parallel, and algebraic
multigrid method, easily adjusted to the underlying problem and computer
architecture, even by the non-expert.
It is very interesting to incorporate information available from the approximate
inverses into the coarsening strategy and grid transfer operators, as suggested
in [19,20]. The expected benet would include improved robustness and local
adaptivity for these multigrid components as well. The authors are currently
pursuing these issues and will report on them elsewhere in the near future.

Acknowledgment

We thank Klaus Stuben for useful comments and suggestions.



--R

An MPI implementation of the SPAI preconditioner on the T3E
A block version of the SPAI preconditioner
Iterative solution of large sparse linear systems arising in certain multidimensional approximation problems
Frequency domain behavior of a set of parallel multigrid smoothing operators
A sparse approximate inverse preconditioner for the conjugate gradient method
A comparative study of sparse approximate inverse preconditioners


Approximate inverse preconditioners via sparse-sparse iterations
A priori sparsity patterns for parallel sparse approximate inverse preconditioners
Robustness and scalability of algebraic multigrid

Parallel preconditioning with sparse approximate inverses

Iterative Solution of Large Sparse Systems of Equations
Approximate sparsity patterns for the inverse of a matrix and preconditioning
Factorized sparse approximate inverse preconditionings: I.

Numerical experiments with algebraic multilevel preconditioners
A multilevel AINV preconditioner
Parallel adaptive multigrid


Toward an e
Sparse approximate inverse smoother for multi- grid
Introduction to algebraic multigrid
An Introduction to Multigrid Methods
On the robustness of ILU-smoothing
Matrix prolongations and restrictions in a black-box multigrid solver
--TR
On the robustness of Ilu smoothing
Matrix-dependent prolongations and restrictions in a blackbox multigrid solver
Multigrid methods on parallel computersMYAMPERSANDmdash;a survey of recent developments
Factorized sparse approximate inverse preconditionings I
A Sparse Approximate Inverse Preconditioner for the Conjugate Gradient Method
Parallel Preconditioning with Sparse Approximate Inverses
Approximate Inverse Preconditioners via Sparse-Sparse Iterations
Approximate sparsity patterns for the inverse of a matrix and preconditioning
A comparative study of sparse approximate inverse preconditioners
Toward an Effective Sparse Approximate Inverse Preconditioner
A Priori Sparsity Patterns for Parallel Sparse Approximate Inverse Preconditioners
Robustness and Scalability of Algebraic Multigrid
Sparse Approximate Inverse Smoother for Multigrid
Robust Parallel Smoothing for Multigrid Via Sparse Approximate Inverses
Coarse-Grid Selection for Parallel Algebraic Multigrid

--CTR
Michele Benzi, Preconditioning techniques for large linear systems: a survey, Journal of Computational Physics, v.182 n.2, p.418-477, November 2002
robust smoothing;algebraic multigrid;approximate inverses;parallel multigrid