Robustness and Scalability of Algebraic Multigrid.
Algebraic multigrid (AMG) is currently undergoing a resurgence in popularity, due in part to the dramatic increase in the need to solve physical problems posed on very large, unstructured grids. While AMG has proved its usefulness on various problem types, it is not commonly understood how wide a range of applicability the method has.  In this study, we demonstrate that range of applicability, while describing some of the recent advances in AMG technology.  Moreover, in light of the imperatives of modern computer environments, we also examine AMG in terms of algorithmic scalability. Finally, we show some of the situations in which standard AMG does not work well and indicate the current directions taken by AMG researchers to alleviate these difficulties.
Introduction
. Algebraic multigrid (AMG) was rst introduced in the early
1980's [11, 8, 10, 12], and immediately attracted substantial interest [32, 28, 30, 29].
Research continued at a modest pace into the late 1980's and early 1990's [18, 14, 21,
25, 20, 26, 22]. Recently, however, there has been a major resurgence of interest in the
eld, for \classical" AMG as dened in [29], as well as for a host of other algebraic-
type multilevel methods [3, 16, 34, 6, 2, 4, 5, 15, 33, 17, 35, 36, 37]. Largely, this
resurgence in AMG research is due to the need to solve increasingly larger systems,
with hundreds of millions or billions of unknowns, on unstructured grids. The size
of these problems dictates the use of large-scale parallel processing, which in turn
demands algorithms that scale well as problem size increases. Two dierent types of
scalability are important. Implementation scalability requires that a single iteration
be scalable on a parallel computer. Less commonly discussed is algorithmic scalability,
which requires that the computational work per iteration be a linear function of the
problem size and that the convergence factor per iteration be bounded below 1 with
bound independent of problem size. This type of scalability is a property of the
algorithm, independent of parallelism, but is a necessary condition before a scalable
implementation can be attained.
Multigrid methods are well known to be scalable (both types) for elliptic problems
on regular grids. However, many modern problems involve extremely complex
geometries, making structured geometric grids extremely di-cult, if not impossible,
to use. Application code designers are turning in increasing numbers to very large
unstructured grids, and AMG is seen by many as one of the most promising methods
for solving the large-scale problems that arise in this context.
This study has four components. First, we examine the performance of \classical"
AMG on a variety of problems having regular structure, with the intent of determining
its robustness. Second, we examine the performance of AMG on the same suite of
problems, but now with unstructured grids and/or irregular domains. Third, we
study the algorithmic scalability of AMG by examining its performance on several of
Center for Applied Scientic Computing (CASC), Lawrence Livermore National Laboratory,
Livermore, CA. Email:fcleary, rfalgout, vhenson, jjonesg@llnl.gov
y Department of Applied Mathematics, University of Colorado, Boulder, CO. Email: ftmanteuf,
stevemg@boulder.colorado.edu
z USS Florida (SSBN-728), Naval Submarine Base, Silverdale, WA, Email: JerryTrish@aol.com
x Front Range Scientic, Boulder, CO. Email: jruge@sobolev.Colorado.EDU
the problems using grids of increasing sizes. Finally, we introduce a new method for
computing interpolation weights, and we show that in certain troublesome cases it
can signicantly improve AMG performance.
Our study diers from previous reports on the performance of AMG (e.g., [29, 30])
primarily by our examination of algorithmic scalability, our emphasis on unstructured
grids, and the introduction of a new algorithm for computing interpolation weights.
In Section 2, a description of some details of the AMG algorithm is given to provide
an understanding of the results and later discussion. In Section 3, we present results
of AMG applied to a range of symmetric scalar problems, using nite element discretizations
on structured and unstructured 2D and 3D meshes. AMG is also tested
on nonsymmetric problems, on both structured and unstructured meshes, and the results
are presented in Section 4. A version of AMG designed for systems of equations
is tested, with the focus on problems in elasticity. Results are discussed in Section
5. In Section 6, we introduce and report on tests of a new method for computing
interpolation weights. We concluding with some remarks in Section 7.
2. The Scalar AMG Algorithm. We begin by outlining the basic principles
and techniques that comprise AMG. Detailed explanations may be found in [29].
Consider a problem of the form
(1)
where A is an n  n matrix with entries a ij . For convenience, the indices are iden-
tied with grid points, so that u i denotes the value of u at point i, and the grid is
denoted
by
ng. In any multigrid method, the central idea is that error
e not eliminated by relaxation must be removed by coarse-grid correction. Applied to
elliptic problems, for example, simple relaxations (Jacobi, Gauss-Seidel) reduce high
frequency error components e-ciently, but are very slow at removing smooth compo-
nents. However, the smooth error that remains after relaxation can be approximated
accurately on a coarser grid. This is done by solving the residual equation
on a coarser grid, then interpolating the error back to the ne grid and using it to
correct the ne-grid approximation. The coarse-grid problem itself is solved by a recursive
application of this method. One iteration of this process, proceeding through
all levels, is known as a multigrid cycle. In geometric multigrid, standard uniform
coarsening and linear interpolation are often used, so the main design task is to choose
a relaxation scheme that reduces errors the coarsening process cannot approximate.
One purpose of AMG is to free the solver from dependence on geometry, so AMG
instead xes relaxation (normally Gauss-Seidel), and its main task is to determine a
coarsening process that approximates error that this relaxation cannot reduce.
An underlying assumption in AMG is that smooth error is characterized by small
residuals, that is, Ae  0, which is the basis for choosing coarse grids and dening
interpolation weights. For simplicity of discussion here, we assume that A is a symmetric
positive-denite M-matrix, with a ii > 0; a ij  0 for j 6= i, and
This assumption is made for convenience; AMG will frequently work well on matrices
that are not M-matrices. To dene any multigrid method, several components are
required. Using superscripts to indicate level number, where 1 denotes the nest level
so that A
and
the components that AMG needs are as follows:
1.
\Grids"
M .
2. Grid operators A 1 ; A
3. Grid transfer operators:
Interpolation I k
Restriction I k+1
4. Relaxation scheme for each level.
Once these components are dened, the recursively dened cycle is as follows:
Algorithm: MV k
Otherwise:
times on A k u
Perform coarse grid correction:
Set
\Solve" on level k+1 with MV k+1
Correct the solution by
times on A k u
For this cycle to work e-ciently, relaxation and coarse-grid correction must work
together to eectively reduce all error components. This gives two principles that
guide the choice of the components:
P1: Error components not e-ciently reduced by relaxation must be
well approximated by the range of interpolation.
P2: The coarse-grid problem must provide a good approximation to
ne-grid error in the range of interpolation.
Each of these aects a dierent set of components: given a relaxation scheme,
P1 determines the coarse grids and interpolation, while P2 aects restriction and
the coarse grid operators. In order to satisfy P1, AMG takes an algebraic
relaxation is xed, and the coarse grid and interpolation are automatically chosen so
that the range of the interpolation operator accurately approximates slowly diminishing
error components (which may not always appear to be \smooth" in the usual
sense). P2 is satised by dening restriction and the coarse-grid operator by the
Galerkin formulation:
I k+1
and A
(2)
When A is symmetric positive denite, this ensures that the correction from the
exact solution of the coarse-grid problem is the best approximation in the range of
interpolation [23], where \best" is meant in the A-norm: by jjvjj A  hAv; vi 1=2 .
The choice of components in AMG is done in a separate preprocessing step:
1.
2.
Partition
into disjoint sets C k and F k .
(a)
(b) Dene interpolation I k
.
3. Set I k+1
and A
.
4.
If
is small enough, set
set go to step 2.
Step 2 is the core of the AMG setup process. Since the focus is on coarsening a
particular level k, such superscripts are omitted here and c and f are substituted for
necessary to avoid confusion. The goal of the setup phase is to
choose the set C of coarse-grid points and, for each ne-grid point
small set C i  C of interpolating points. Interpolation is then of the form:
I f
c u c
2.1. Dening Interpolation Weights. To dene the interpolation weights
recall that slow convergence is equivalent to small residuals, Ae  0. Thus,
we focus on errors satisfying
a ii e i
Now, for any a ij that is relatively small, we could substitute e i for e j in (4) and
this approximate relation would still hold. This motivates the denition of the set
of dependencies of a point i, denoted by S i , which consists of the set of points j for
which a ij is large in some sense. Hence, i depends on such j because, to satisfy the ith
equation, the value of u i is aected more by the value of u j than by other variables.
The denition used in AMG is
( a ik )
with  typically set to be 0.25. We also dene the set S T
that is, the
set of points j that depend on point i, and we say that
i is the set of in
uences of
point i. Note: our terminology here diers from the classical use in [29], which refers
to i as being strongly connected to or strongly dependent on j if j 2 S i and which uses
no specic terminology for
A basic premise of AMG is that relaxation smoothes the error in the direction of
in
uence. Hence, we may select C as the set of interpolation points for i,
and adhere to the following criterion while choosing C and F :
P3: For each i 2 F , each j 2 S i is either in C or S j \ C i 6= ;.
That is, if i is a ne point, then the points in
uencing i must either be coarse points
or must themselves depend on the coarse points used to interpolate u i . This allows
approximations necessary to dene interpolation. For can be rewritten as:
a ii e i
a ik e k
AMG interpolation is dened by making the following approximation in (6):
a jk e k
a jk
otherwise.
Substituting this into (6) and solving for e i gives the desired interpolation weights for
point
2.2. Selecting the Coarse Grid. The coarse grid is chosen to satisfy the criterion
above, while attempting to control its size. We employ the two-stage process
described in [29], modied slightly to re
ect our modied terminology. The grid is
rst \colored", providing a tentative C/F choice. Essentially, a point with the largest
number of in
uences (\in
uence count") is colored as a C point. The points depending
on this C point are colored as F points. Other points in
uencing these F points
are more likely to be useful as C points, so their in
uence count is increased. The
process is repeated until all points are either C or F points.
Details of the initial C=F choice are as follows:
Repeat until
k . Set
(the number of
points depending on the point i).
Pick an i 2 U with maximal  i . Set
T fig and
fig.
For
(points depending on fig) do:
S fjg and
For all k 2 S j
T U set  (Increment the
for points that in
uence the new F -points).
For
T U set
Next, a second pass is made, in which some F points may be recolored as C points
to ensure that P3 is satised. In this pass, each F -point i is examined. The coarse
interpolatory set C
S C is dened. Then, if i depends on another F -point, j,
the points in
uencing j are scanned, to see if any of them are in C i . If this is not the
case then j is tentatively converted into a C-point and added to C i . The dependencies
of i are then examined anew. If all F -points depending on i now depend on a point
in C i then j is permanently made a C-point and the algorithm proceeds to the next
F -point and repeats. If, however, the algorithm nds another F -point dependent on
i that is not dependent on a point in C i then i itself is made into a C-point and j
returned to the pool of F -points. This procedure is followed to minimize the number
of F -points that are converted into C-points.
We make a brief comment about the computational and storage costs of the
setup phase. Unlike geometric multigrid, these costs cannot be predicted precisely.
Instead, computational cost must be estimated based on the average \stencil size"
over all grids, the average number of interpolation points per F -point, the ratio of the
total number of gridpoints on all grids to the number of points on the ne grid (grid
complexity), and the ratio of the number of nonzero entries in all matrices to that of
the ne-grid matrix (operator complexity). While a detailed analysis is beyond the
scope of this work [29], a good rule of thumb is that the computational eort for the
setup phase is typically equivalent to between four and ten V -cycles.
3. Results for Symmetric Problems. In this section, results for AMG applied
to symmetric scalar problems are presented. Initially, constant-coe-cient diusion
problems in 2D are tested as a baseline for comparison as we begin to introduce
complications, including unstructured meshes, irregular domains, and anisotropic and
discontinuous coe-cients. Results for 3D problems follow. All problems are run using
the same AMG solver with xed parameters. On many problems, it is possible to
improve our results by tuning some of the input parameters (there are many), but
the purpose here is to show AMG's basic behavior and robustness over a range of
problems.
The primary indicator of the speed of the algorithm is the asymptotic convergence
factor per cycle. This is determined by applying 20 cycles to the homogeneous prob-
lem, starting with a random initial guess, then measuring the reduction in the norm
of the residual from one cycle to the next (we use the homogeneous problem to avoid
contamination by machine representation). Generally, this ratio starts out very small
for the rst few cycles, then increases to some asymptotic value after 5-10 cycles,
when the most slowly converging components become dominant. This asymptotic
value is also a good indicator of the actual error reduction from one cycle to the next.
We use the 2-norm of the residual, although it is easy to show that the asymptotic
convergence factor is just the spectral radius of the AMG V -cycle iteration operator,
and hence is independent of the choice of norm.
The times given are for the setup and a single (1,1) V -cycle. Setup time is what
it takes to choose the coarser grids, dene interpolation, and compute the coarse grid
matrices. Cycle time is for one cycle, not the full solution time. Three machines are
used in this study. The majority of the smaller tests are performed on a Pentium
166MHz PC, although some are performed on a Sun Sparc Ultra 1. For the larger
problems that demonstrate scalability, we use a DEC Alpha. For this reason, timings
should be compared only within individual problems. Additionally, timings for the
smallest problems can have a high relative error, so the larger tests should give a
better picture of performance. Grid complexity is dened as
the number of grid points on level k. This gives an idea of how quickly the grids
are reduced in size. For comparison, in standard multigrid, the number of points is
reduced by a factor of 4 in 2D and 8 in 3D, yielding grid complexities of 4/3 and
8/7, respectively. AMG tends to coarsen more slowly. Operator complexity, which is
a better indicator of the work per cycle, is dened as
is the
average number of non-zero entries per row (or \stencil size") on level k. Thus, the
operator complexity is the ratio of the total number of nonzero matrix entries on all
levels to those on the nest level. Since relaxation work is proportional to the number
of matrix entries, this gives a good idea of the total amount of work in relaxation
relative to relaxation work on the nest grid, and also of the total storage needed
relative to that required for the ne grid matrix. In geometric multigrid, the grid and
operator complexities are equal, but in AMG, operator complexity is usually higher
since average stencil sizes tend to grow somewhat on coarser levels. Note that the
convergence factors and complexities are entirely independent of the specic machine
on which a test is performed.
In the tests reported here, the focus is on nite element discretizations of
r
d 21 (x; y) d 22 (x; y)
Several dierent meshes and diusion coe-cients D are used.
3.1. Regular domains, structured and unstructured grids. The rst ve
problems are 2D Poisson equations, with d
Dierent domains and meshes are used to demonstrate the behavior of AMG with
simple equations.
We begin with the simplest 2D model problem. The success of AMG on the
regular-grid Poisson problem is well-documented [30, 28, 29], so our purpose here is
more to assess its scalability.
Problem 1 This is a simple 5-point Laplacian operator with homogeneous Dirichlet
boundary conditions on the unit square. The experiment is run for uniform meshes
Convergence factor vs. log of problem size
Cycle Times vs. log of problem size
vs. log of problem size
Fig. 1. Top Left: Convergence factors, as a function of number of mesh points, for Problem 1
the uniform-mesh 5-point Laplacian. Top Right: Log-log plots of setup times (circles) and cycle times
(triangles) for the uniform-mesh 5-point Laplacian. The dotted line, for reference, shows perfectly
linear scaling. Bottom: Operator (circles) and grid (triangle) complexities for the uniform-mesh
5-point Laplacian.
with n  n interior grid points, yielding mesh sizes
90000, 250000, and 490000.
Results for Problem 1 are displayed in Figure 1. The convergence factor (per cy-
cle) is very stable at approximately 0.04 for all problem sizes. Both the setup and cycle
time are very nearly linear in N (compare with the dotted line depicting a perfectly
linear hypothetical data set). Here, setup time averages roughly the time of 6 cycles.
As noted before, the operator complexities are higher than the corresponding grid
complexities, but both appear to be unaected by problem size. These data indicate
that AMG (applied to the uniform-mesh Laplacian) is algorithmically scalable: the
computational work is O(N) per cycle and the convergence factor is O(1) per cycle.
An important component of our study is to determine to what extent this algorithmic
scalability is retained as we increase problem complexity.
Problem 2 This is the same equation as Problem 1 ( discretized
on an unstructured triangular mesh. These meshes are obtained from uniform triangulations
by randomly choosing 15-20% of the nodes and \collapsing" them to neighboring
nodes, then smoothing the resulting mesh. The resulting operators might be
represented by M-matrices in some cases, but this is not generally the case. We use
meshes with typical example is shown at
top left in Figure 2.
Results of the experiments are displayed in Figure 2. On the unstructured meshes,
convergence factors tend to show some dependence on mesh size, growing to around
0.35 on the nest grid. It should be noted, however, that these grids tend to be less
structured than many found in practice, and no care was taken to ensure a \good"
mesh; the meshes may have diering characteristics (such as aspect ratios), as there is
a large degree of randomness in their construction. Complexities are also higher with
the unstructured meshes, and the setup time increases correspondingly. The main
point here is that AMG can deal eectively with unstructured meshes without too
much degradation in convergence over the uniform case.
3.2. Irregular domains. We continue to use the Laplacian, but now with irregular
domains. Since our emphasis here is the eects of this irregularity, we restrict our
tests to two representative mesh sizes that give just a snapshot of algorithm scalability.
Problem 3 The computational domain is an unstructured triangular discretization
of the torus 0:05
dierent mesh sizes were used, resulting in
grids with conditions around the hole are
imposed, with Neumann conditions on the outer boundary.
Problem 4 The domain for this problem is shown in Figure 3. The boundary conditions
are Neumann except that a Dirichlet condition is imposed around the small
hole on the right. The meshes are uniform, with resulting in
meshes with respectively. The domain does not easily admit
much coarser meshes.
Problem 5 The domain for this problem is shown on the bottom in Figure 3. Dirichlet
conditions are imposed on the exterior boundary, and Neumann conditions are on
the interior boundaries. A triangular unstructured mesh is used.
Results for Problems 3{5 are given in Table 1. Among these problems, Problem 3
has the simplest domain, but the least structured mesh and the slowest convergence.
This indicates that domain conguration generally has little eect on AMG behavior,
while the structure (and perhaps the quality) of the mesh is more important.

Table

1. Results for Problems 3{5.
Poisson problem on unstructured meshes, irregular domains
Convergence
Problem N factor/cycle (sec) (sec) complexity complexity
Convergence Factor vs. log of problem size
Cycle Times vs. log of problem size
vs. log of problem size
Fig. 2. Top Left: A typical unstructured grid for Problem 2, obtained by randomly deleting
15% of the nodes in a regular grid and smoothing the result. Top Right: Convergence factors, as
a function of number of mesh points, for the unstructured-mesh 5-point Laplacian. Bottom Left:
Log-log plots of setup times (circles) and cycle times (triangles) for the unstructured-grid 5-point
Laplacian. The dotted line, for reference, shows perfectly linear scaling. Bottom Right: Operator
(circles) and grid (triangle) complexities for the unstructured-grid 5-point Laplacian.
Fig. 3. Domain (Top Left) and typical grid (Top Right) for Problem 4. Note that the mesh
size necessary to display the triangulation is too coarse to observe the Dirichlet hole. Finer meshes
are used for the calculations. Bottom: Typical grid for Problem 5.
3.3. Isotropic diusion. The next problem set deals with isotropic diusion:
Discontinuous d(x; y) can cause problems for
many solution methods, including standard multigrid methods, although it is possible
to get good results either by aligning the discontinuities along coarse grid lines, or
by using operator-dependent interpolation [1]. In AMG, nothing special is required,
since it is based on operator-dependent interpolation. The problems are categorized
according to the diusion coe-cient used. The unit square is discretized on four
meshes: two structured meshes with and two unstructured
meshes, with 54518. The diusion coe-cients are dened in terms of
a parameter, c, allowed to be either 10 or 1000, as follows:
Problem 6 d(x;
Problem 7 d(x;
Problem 8 d(x;
1:0 0:125  max (jx 0:5j; jy 0:5j)  0:25;
c otherwise:
Problem 9 d(x;
1:0 0:125
c otherwise:
Results for these problems are presented in Table 2, which contains observed
convergence factors and operator complexities for the various combinations of grid
size and type, diusion coe-cient function, and discontinuity jump size. The overall
results are fairly predictable. Convergence factors are fairly uniform. On the structured
meshes, they tend to grow slightly with increasing grid size. They are noticeably
larger for unstructured grids, and they appear to grow somewhat with increasing grid
size. (As noted before, comparison among unstructured grids of various sizes must
take into account that their generation involves some randomness, so they may dier
in important ways.) The convergence factor does not seem to depend signicantly
on the size of the jump in the diusion coe-cient. In many cases, results were better
with Indeed, AMG has been applied successfully
to problems with much larger jumps [28]; see also Problem 17. Note that there are
only minor variations in operator complexity for the dierent problems and dierent
grid sizes. The only signicant eect on operator complexity appears to be whether
the grid is structured or unstructured, with the latter showing complexity increases
of about 30{40%. It should be noted, however, that even in these cases, the entire
operator hierarchy can be stored in just over three times the storage required for the
ne-grid alone.
Table

2. Results for Problems 6{9.
Poisson problem, variable and discontinuous coe-cients
Uniform mesh size Unstructured mesh size
Problem 16642 66049 13755 54518
# c conv. cmplxty conv. cmplxty conv. cmplxty conv. cmplxty
6 1000 0.097 2.21 0.180 2.2 0.264 3.32 0.369 3.36
9 1000 0.171 2.35 0.168 2.30 0.234 3.31 0.298 3.40
Problems 10{13 are designed to examine the case in which the diusion coe-cient
is discontinuous and to determine whether the \scale" of the discontinuous regions
aects performance. Accordingly, Problems 10{12 use a \checkerboard" pattern:
even and i
where Specically,
Problem 2.
Problem
Problem
For the last problem of this group, we have
Problem 13 d(x;
Results for Problems 10{13 are displayed in Table 3. The overall trend is similar
to the results for Problems 6{9, showing convergence factors that grow slightly with
problem size and that are noticeably larger for unstructured grids.

Table

3. Results for Problems 10{13.
Poisson problem, variable and discontinuous coe-cients
Uniform mesh size Unstructured mesh size
Problem 16642 66049 13755 54518
# c conv. cmplxty conv. cmplxty conv. cmplxty conv. cmplxty
The best results for the isotropic diusion problems are obtained for Problem 6,
where the coe-cient is continuous. The worst convergence factor is obtained for the
50  50 \checkerboard" pattern of Problem 12. Note that AMG performs well on
Problem 8, where \smooth" functions are approximately constant in the center of the
region, zero in the high-diusion zone near the boundary, and smoothly varying in be-
tween. Good interpolation in the low-diusion band is essential to good convergence.
AMG appears to work quite well with discontinuous diusion coe-cients,
even when they vary randomly by a large factor from point to point, as in Problem
13.
3.4. Anisotropic diusion. The next series of problems deals with anisotropic
diusion, which can arise in several ways. Anisotropy can be introduced by the mesh
being rened dierently in each directions, perhaps to resolve a boundary layer or
some other local phenomenon. Another case is a tensor product grid used in order to
rene some area [x with the mesh size small for x 2 [x
large elsewhere. This maintains a logically rectangular mesh, but causes
anisotropic discretizations in dierent parts of the domain. This is relatively easy to
deal with in geometric multigrid, where line relaxation and/or semi-coarsening can
be used [9]. Non-aligned anisotropy, which is more di-cult to handle with standard
multigrid, arises from the operator itself, such as the case of the full potential operator
in transonic
ows. The performance of AMG on grid-induced (aligned) anisotropy
has been reported previously [29], so we instead focus here on non-aligned anisotropy.
Both types of anisotropy can be written in terms of the diusion equation using the
coe-cient matrix:
cos  sin  sin 2
When  is constant, this gives the operator @  +@  , where  is in the direction . On
a rectangular grid with mesh sizes h x and h y , the usual Poisson equation corresponds
to the diusion equation with
x .
Problem 14 This problem features a non-aligned anisotropic operator on the unit
square, with Dirichlet boundary conditions at conditions
on the other two sides. The cases were both examined
with and =4. Each such combination is discretized on a uniform
square mesh (with bilinear elements) and both uniform and unstructured triangular
meshes (with linear elements). The uniform meshes have
the unstructured meshes have
Convergence factors for Problem 14, as shown in Figure 4, generally degrade with
increasing . This is to be expected, as it indicates lessened alignment of anisotropy
with the grid directions. The strong anisotropy case yields convergence factors as
high as 0.745. As noted above, the non-aligned case is very di-cult, even for standard
multigrid, and is the subject of ongoing study. One encouraging result is that
the unstructured grid formulations are relatively insensitive to grid anisotropy, with
convergence factors that hover between 0.3 and 0.5 in all cases except
these results indicate that AMG is rather robust for anisotropic problems, although
convergence factors are somewhat higher than those typically obtained with AMG on
isotropic problems.
Convergence factor vs.  Convergence factor vs.
0.80.20.6Fig. 4. Convergence factors, plotted as a function of anisotropy direction . Left: Moderate
anisotropy, In each plot, the solid lines are the
uniform meshes with square discretizations, the dotted lines are uniform meshes with triangular
discretizations, and the dashed lines are the unstructured triangulations. In each case, the larger
grid size is indicated by a symbol (\o", \+", or \").
Problem 15 We use the operator of Problem 4, r  (Dru), but with
This yields a discretization such that, on any circle
centered at the origin, there are dependencies in the tangential direction, but none in
the radial direction.
This problem is very di-cult to solve by conventional methods. Using the same
meshes as in Problem 14, AMG produced the convergence factors given in Table 4.

Table

4. Results for Problem 15: Circular diusion coe-cient
uniform mesh (square) uniform mesh (triangular) unstructured mesh
The convergence factors illustrate the di-culty with this problem, which cannot be
handled easily by geometric methods, even on regular meshes. A polar-coordinate
mesh would allow block relaxation over strongly coupled points, but would suer from
the di-culties of polar-coordinate grids (e.g., singularity at the origin) and would be
useless for more general anisotropies. While convergence of AMG here is much slower
than what we normally associate with multigrid methods, this example shows that
AMG can be useful even for extremely di-cult problems.
3.5. 3D problems. Turning our attention to three dimensions, we do not expect
special di-culties here, since AMG is based on the algebraic relationships between
the variables.
Problem This is a 3D Poisson problem on the unit cube. Discretization is by
trilinear nite elements on a rectangular mesh. Dirichlet boundary conditions are
imposed at conditions are imposed at the other
boundaries. Mesh line spacing and the number of mesh intervals were both varied
to produce several grids with dierent spacings and extents in the three coordinate
directions. The various combinations of mesh sizes used, convergence factors, and
operator complexities are shown in Table 5.

Table

5. Results for Problem 16.
3D Poisson problem, regular rectangular mesh
Convergence Operator
Nx hx Ny hy Nz hz N factor/cycle Complexity
Problem 17 This is a 3D unstructured mesh problem, generated by a code used at
Lawrence Livermore National Laboratory, for the diusion problem
g(~x). The domain is a segment of a sphere, from
and to =2. The coe-cient a(~x) is a large constant for r  0:05 and a small
constant for r > 0:05, with a step discontinuity of 1:010 26 . The boundaries
and are Dirichlet, while the boundaries are surfaces of
symmetry. Discretization is by nite elements using hexahedral elements.

Figure

5 shows the locations of the nodes at the element corners for one of the prob-
lems. Three problem sizes are given, with Convergence
factors and operator complexities for these problems are given in Table 6.

Table

6. Results for Problem 17.
3D unstructured diusion problem
Convergence Operator
8000 0.166 2.84
AMG apparently works quite well for 3D problems, including those with discontinuous
coe-cients. The convergence factors are good in all cases. Complexity varies
signicantly, with the highest values for the uniform grids, but decreasing markedly
with increasing grid anisotropy. This may be taken as further evidence that AMG
automatically takes advantage of directions of in
uence.
performed well on this suite of symmetric scalar test problems.
Many of these problems are designed to be very di-cult, often unrealistically so,
especially those with the circular anisotropic diusion pattern and the random diu-
sion coe-cients. Recall that the same AMG algorithm, with no parameter tuning,
was used in all cases. There are a number of tools for increasing the e-ciency of
AMG, especially on symmetric problems, that have proved useful in many cases. One
0.10.010.030.050.07Fig. 5. Nodes at element corners, 3D diusion problem. Top Left: View of all the node
locations. Top Right: View from directly overhead, showing radial lines of nodes in the azimuthal
direction. Bottom: Distribution of nodes within a plane of constant azimuth.
is the so-called V ? cycle [30], in which the coarse-grid corrections are multiplied by
an optimal parameter, determined by minimizing the A-norm of the corrected error.
Another, which has been successful in applying AMG to Maxwell's equations [27], is
to use an outer conjugate gradient iteration, with AMG cycling as a preconditioner.
Often, when AMG fails to perform well, the problem lies in a small number of components
that are not reduced e-ciently by relaxation or coarse-grid correction, and
conjugate gradients can be very e-cient in such cases. Other methods for improving
e-ciency include the F cycle [7, 31] and the full multigrid (FMG) method, whose
applicability to AMG is the subject for future research.
4. AMG applied to Nonsymmetric Scalar Problems. Although much of
the motivation and theory for AMG is based on symmetry of the matrix, this is not
at all a requirement for good convergence behavior. Mildly nonsymmetric problems
behave essentially like their symmetric counterparts. Such cases arise when a nonsymmetric
discretization of a symmetric problem is used or when the original problem is
predominantly elliptic. An important requirement for current versions of AMG is that
point Gauss-Seidel relaxation converge, however slowly. Thus, central dierencing of
rst-order terms, when they dominate, cannot be used because of severe loss of diagonal
dominance. Even in these cases, successful versions of AMG can be developed
using Kaczmarz relaxation [9]. Nevertheless, we restrict ourselves here to upstream
dierencing so that we can retain our use of Gauss-Seidel relaxation.
Problem This is a convection-diusion problem of the form
sin
with Dirichlet boundary conditions. Triangular meshes are used, both structured
unstructured 54518). The diusion term
is discretized by nite elements. The convection term is discretized using upstream
dierencing, that is, the integral of the convection term is computed over the triangle
and added to the equation corresponding to the node with the largest coe-cient (the
node \most upstream"). Note that this can result in a matrix that has o-diagonal
entries of both signs. Two choices for  are employed:
are conducted with
Results are presented in Figure 6. The curves in the top left graph are for the
structured grid results are displayed with solid lines, and the unstructured-grid results
are displayed with dashed lines. For each pair of curves, the curve with the marker (\o"
or \") indicates the mesh with larger N . The convection-dominated case
is shown at the top right (structured grids) and on the bottom (unstructured grids).
In each case, the smaller N is shown with solid lines and the larger N with dashed
lines. Note that convergence is generally good and fairly uniform, particularly for the
unstructured cases. Results on the smaller uniform mesh are especially good when the
ow is aligned with the directions or =2. This is due to the triangulation:
to obtain the uniform mesh, the domain is partitioned into squares, and then each
square is split into two triangles, with the diagonal going from the lower left to the
upper right; the \good" directions are aligned with the edges of the triangles. This
also has an eect on the quality of the discretization, and on convergence, when the
ow is in the direction 3=4 and 7=4. Here, the discretization used for the convection
term causes a rather severe loss of positivity in the o-diagonals. This is more the
fault of the discretization than AMG. For with the smaller uniform mesh,
Convergence factor vs.  Convergence factor vs.
structured grid
Convergence factor vs.
unstructured grid
Fig. 6. Convergence factors, plotted as a function of direction, , for the nonsymmetric prob-
lem. Top Left: Weaker convection case, 0:1 The solid lines are the uniform meshes, and the
dashed lines are the unstructured triangulations. In each case the larger grid size is indicated by
the placement of a symbol (\o" or \"). Top Right: Convection-dominated case, structured
grid. The solid line is the smaller N , the larger N is indicated by the dashed line. Bottom:
Convection-dominated case, unstructured grid. The solid line is the smaller N , the
larger N is indicated by the dashed line.
AMG was unable to handle the discretization produced and failed in the setup phase.
Often in multigrid applications, problems in convergence indicate problems with the
discretization, as is the case here. Note that for unstructured grids where
ow cannot
align with (or against) the grid, convergence is generally more uniform. An interesting
point is that, in many cases, a smaller diusion coe-cient reduces the convergence
factor. This is particularly striking with the unstructured meshes. Finally, note that
there is generally not much dierence between results for  and for  + , so there
is no benet from accidental alignment of relaxation with the
ow direction, and,
conversely, there is no slowing of convergence due to upstream relaxation. This is
due to the C/F ordering of relaxation. These tests show that AMG can be applied
to nonsymmetric problems. While the convergence factors in this test are generally
less than 0:2{0:25, which is certainly acceptable, some concern may be raised about
the scalability issue, since the convergence factors for the larger unstructured mesh
are noticeably greater than those for the smaller unstructured grids. It remains to be
determined if the convergence factors continue to grow with increasing problem size,
or if they reach an asymptotic limit.
5. AMG for Systems of Equations. The extension of AMG to \systems"
problems, where more than one function is being approximated, is not straightforward.
Many dierent approaches can be formulated. Consider a problem with two unknown
functions of the form
u
The scalar algorithm could work in special circumstances (for example, if B and C are
relatively small in some sense), but, generally, the scalar ideas of smoothness break
down. One approach would be to iterate in a block fashion on the two equations,
with two separate applications of AMG, one using A as the matrix (solving for u,
holding v xed) and one using B (solving for v, holding u xed), and repeating until
convergence. This is often very slow. An alternative is to use this block iteration as
a preconditioner in an outer conjugate gradient solution process.
Another fairly simple alternative is to couple the block iteration process on all
levels, that is, to coarsen separately for each function, obtaining two interpolation
operators I u and I v , then dene a full interpolation operator of the form
I
The Galerkin approach can then be used to construct the coarse grid operator,
I T
I T
Once the setup process is completed, multigrid cycles are performed as usual. We will
call this the function approach since it treats each function separately in determining
coarsening and interpolation. When u and v are dened on the same grid, it is
also possible to couple the coarse-grid choices for both, allowing for nodal relaxation,
where both unknowns are updated simultaneously at a point. Following are results
for the function approach applied to several problems in 2D and 3D elasticity.
Problem 19 This problem is plane-stress elasticity:
where u and v are displacements in the x and y directions, respectively. This can
be a di-cult problem for standard multigrid methods, especially when the domain is
long and thin. The problem is discretized on a rectangular grid using bilinear nite
elements, and several dierent problem sizes and domain congurations are used.
Problem 20 This problem is 3D elasticity:
where u, v, and w are displacements in the three coordinate directions. The problem
is discretized on a 3D rectangular grid using trilinear nite elements. Several dierent
problem sizes and domain congurations are used.
In all tests, we take 0:3. The function approach with (1,1) V -cycles is used in
all tests. Results for Problem 19 are contained in Table 7, and for Problem 20 in Table
8. Note that complexities are stable in 2D, with some dependence on problem size
in 3D. Convergence depends fairly heavily on the number of xed boundaries in both
2D and 3D, with convergence degrading as the number of free boundaries increases.

Table

7. Results for 2D elasticity, Problem 19
# of xed Convergence Operator

Table

8. Results for 3D elasticity, Problem 20
# of xed Convergence Operator
6. Iterative Interpolation Weights. Occasionally, we encounter situations
where convergence of AMG is poor, yet no specic reason is apparent. Our experience
leads us to believe that the fundamental problem, in many cases, stems from the
limitation of the matrix entry a jk to re
ect the true "smoothness" between e j and
e k . Often the true in
uences between variables are not clear. One case where this
limitation is quite evident is where nite elements with extreme aspect ratios are
used, especially in cases of extreme grid anisotropy or of thin-body elasticity. As a
simple example, consider the 2D nine-point negative Laplacian based on quadrilateral
elements that are stretched in the x-direction. The stencil changes as
y !1:
The limiting case is no longer an M-matrix. Indeed, even moderate aspect ratios
(e.g., o-diagonal entries of both signs. It is not immediately clear
how the neighbors to the east and west of the central point should be treated. Do
they in
uence the central point? Should they be in S i ? Even if they are not treated
as in
uences, similar questions arise about how the corner points relate to the central
point. Geometric intuition indicates they are decoupled from the center, and should
not be treated as in
uences. Yet, for the most common choice of  in (5), AMG treats
them as in
uences. Another di-culty arises when two F points i and j in
uence each
other. Then e j must be approximated in the second sum on the right side of (6) to
determine the weights for i, while e i must be approximated, in (6) but with the roles
of i and j reversed, to determine the weights for j. However, since both e i and e j are
to be interpolated (being F -point values), it makes sense to use the interpolations to
obtain these approximations, that is, the approximations for e i and e j in (6) should
be
and e i
respectively. Note that the approximations for any points in C are unchanged in these
equations.
This gives an implicit system for the interpolation weights, which is solved by
an iterative scheme with the initial approximation . The new interpolation
weights are then calculated in a Gauss-Seidel-like manner, using the most recently
computed weights to make the approximations in (12). Two sweeps are generally
su-cient. An important addition to the process is that, after the rst sweep, the
interpolation sets are modied by removing from C i any point for which a negative
interpolation weight is computed. The second sweep is then used to compute the nal
interpolation weights.
We present results of two experiments illustrating the eectiveness, on certain
types of problems, of using this iterative weight denition scheme. Other examples
may be found in [24].
Problem 21 This operator is the \stretched quadrilateral" Laplacian mentioned
above, discretized on an 10000. The
stretching factor " represents the ratio of the x-dimension to the y-dimension of the
quadrilateral. Values used are In each case, the convergence
factor is computed for the choices 0:5. In the latter case, the corner
points are not treated as in
uences, and AMG selects a semi-coarsened coarse grid,
which is the method geometric multigrid would take.

Table

9. Results for problem 21,
Convergence Factor
standard weights iterative weights
100 900 0.83 0.53 0.82 0.23
100 10000 0.93 0.55 0.93 0.28
Results are displayed in Table 9. On the smallest problem, iterative weights have
no eect. However, the convergence rate on that problem is quite good, even for
standard weights. For moderate stretching (" < 100), the eect of iterative weighting
is to correct for misidentied in
uence (i.e., improvement for and to improve
the results even for correctly identied in
uence 0:5). For extreme stretching,
only the latter eect applies.
Problem 22 Here we use the unstructured 3D diusion operator from Problem 17,
whose grid is displayed in Figure 5. Problem sizes are 8000. The
problem includes a very large jump discontinuity, O(10 26 ), in the diusion coe-cients.

Table

10. Results for problem 22
Convergence Factor
standard weights iterative weights
Results are displayed in Table 10. Again we see that, on the smallest problem, iterative
weights have no eect, but that convergence there is fairly good anyway, even for
standard weights. On the two larger problems, iterative weights produce signicant
improvement for both choices of . Apparently, iterative weighting is countering the
eects of both poor element aspect ratios near the boundaries and jump discontinuities
in identifying in
uences among variables.
On these problems, and similar problems characterized by coe-cient discontinuities
and/or extreme aspect ratios in the elements, iterative weight denition proves
to be quite eective. However, iterative weighting is not always eective at improving
slow AMG convergence, and in a few cases it can actually cause very minor
degradation in performance [24]. We study a new approach in [13], called element
interpolation (AMGe), which has the promise of overcoming the di-culties associated
with poor aspect ratios, misidentied in
uences, and thin-body elasticity, provided
the individual element stiness matrices are available.
7. Conclusions. The need for fast solvers for many types of problems, especially
those discretized on unstructured meshes, is a clear indication that there is a market
for software with the capabilities that AMG oers. Our study here demonstrates the
robustness of AMG as a solver over a wide range of problems. Our tests indicate
that it can be further extended, and that robust, e-cient codes can be developed
for problems that are very di-cult to solve by other techniques. AMG is also shown
to have good scalability on model problems. This scalability does tend to degrade
somewhat with increasing problem complexity, but the convergence factors remain
tractable even in the worst of these situations.



--R

The multi-grid methods for the diusion equation with strongly discontinuous coe-cients
Stabilization of algebraic multilevel iteration
The algebraic multilevel iteration methods - theory and applications
A class of hybrid algebraic multilevel preconditioning methods

Towards algebraic multigrid for elliptic problems of second order
Guide to multigrid development

guide with applications to uid dynamics





Outlines of a modular algebraic multilevel method
Interpolation and related coarsening techniques for the algebraic multigrid method
Additive multilevel-preconditioners based on bi-linear 4interpolation
A note on the vectorization of algebraic multigrid algorithms
Matrix Analysis
Convergence of algebraic multigrid methods for symmetric positive de


Multigrid methods for variational problems: general theory for the V-cycle
Interpolation weights of algebraic multigrid
On smoothing properties of SOR relaxation for algebraic multigrid method
A simple parallel algebraic multigrid
Multigrid methods for solving the time- harmonic Maxwell equations with variable material parameters









An energy-minimizinginterpolation for robust multi-grid methods
--TR

--CTR
Emden Henson , Ulrike Meier Yang, BoomerAMG: a parallel algebraic multigrid solver and preconditioner, Applied Numerical Mathematics, v.41 n.1, p.155-177, April 2002
Chi Shen , Jun Zhang , Kai Wang, Distributed block independent set algorithms and parallel multilevel ILU preconditioners, Journal of Parallel and Distributed Computing, v.65 n.3, p.331-346, March 2005
Randolph E. Bank, Compatible coarsening in the multigraph algorithm, Advances in Engineering Software, v.38 n.5, p.287-294, May, 2007
K. Kraus , J. Schicho, Algebraic Multigrid Based on Computational Molecules, 1: Scalar Elliptic Problems, Computing, v.77 n.1, p.57-75, February 2006
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
Touihri, Variations on algebraic recursive multilevel solvers (ARMS) for the solution of CFD problems, Applied Numerical Mathematics, v.51 n.2-3, p.305-327, November 2004
J. J. Heys , T. A. Manteuffel , S. F. McCormick , L. N. Olson, Algebraic multigrid for higher-order finite elements, Journal of Computational Physics, v.204
Michele Benzi, Preconditioning techniques for large linear systems: a survey, Journal of Computational Physics, v.182 n.2, p.418-477, November 2002
scalability;algebraic multigrid;interpolation;unstructured meshes