An Algebraic Multilevel Multigraph Algorithm.
We describe an algebraic multilevel multigraph algorithm. Many of the multilevel components are generalizations of algorithms originally applied to general sparse Gaussian elimination. Indeed, general sparse Gaussian elimination with minimum degree ordering is a limiting case of our algorithm. Our goal is to develop a procedure which has the robustness and simplicity of use of sparse direct methods, yet offers the opportunity to obtain the optimal or near-optimal complexity typical of classical multigrid methods.
Introduction
. In this work, we develop a multilevel multigraph algorithm.
Algebraic multigrid methods are currently a topic of intense research interest [17, 18,
20, 46, 12, 48, 38, 11, 44, 3, 4, 1, 2, 5, 16, 7, 29, 28, 27, 42, 41, 21]. An excellent recent
survey is given in Wagner [49]. In many \real world" calculations, direct methods are
still widely used [6]. The robustness of direct elimination methods and their simplicity
of use often outweigh the apparent benets of fast iterative solvers. Our goal here is to
try to develop an iterative solver that can compete with sparse Gaussian elimination
in terms of simplicity of use and robustness and to provide the potential of solving
a wide range of linear systems more e-ciently. While we are not yet satised that
our method has achieved this goal, we believe that it is a reasonable rst step. In
particular, the method of general sparse Gaussian elimination with minimum degree
ordering is a point in the parameter space of our method. This implies that in the
worst case, our method defaults to this well-known and widely used method, among
the most computationally e-cient of general sparse direct methods [26]. In the best
case, however, our method can exhibit the near optimal order complexity of the
classical multigrid method.
Our plan is to take well studied, robust, and widely used procedures and data
structures developed for sparse Gaussian elimination, generalize them as necessary,
and use them as the basic components of our multilevel solver. The overall iteration
follows the classical multigrid V-cycle in form, in contrast to the algebraic hierarchical
basis multigraph algorithm developed in [11].
In this work we focus on the class of matrices which are structurally symmetric;
that is, the pattern of nonzeros in the matrix is symmetric, although the numerical
values of the matrix elements may render it nonsymmetric. Such structurally symmetric
matrices arise in the discretizations of partial dierential equations, say, by
the nite element method. For certain problems, the matrices are symmetric and
positive denite, but for others the linear systems are highly nonsymmetric and/or
indenite. Thus in practice this represents a very broad class of behavior. While our
main interest is in scalar elliptic equations, as in the nite element code PLTMG [8],
our algorithms can formally be applied to any structurally symmetric, nonsingular,
sparse matrix.
Sparse direct methods typically have two phases. In the rst (initialization) phase,
Department of Mathematics, University of California at San Diego, La Jolla, CA 92093. The
work of this author was supported by the National Science Foundation under contract DMS-9706090.
y Bell Laboratories, Lucent Technologies, Murray Hill, NJ 07974.
equations are ordered, and symbolic and numerical factorizations are computed. In
the second (solution) phase, the solution of the linear system is computed using the
factorization. Our procedure, as well as other algebraic multilevel methods, also
breaks naturally into two phases. The initialization consists of ordering, incomplete
symbolic and numeric factorizations, and the computation of the transfer matrices
between levels. In the solution phase, the preconditioner computed in the initialization
phase is used to compute solution using the preconditioned composite step conjugate
gradient (CSCG) or the composite step biconjugate gradient (CSBCG) method [9].
Iterative solvers often have tuning parameters and switches which require a certain
level of a priori knowledge or some empirical experimentation to set in any particular
instance. Our solver is not immune to this, although we have tried to keep the number
of such parameters to a minimum. In particular, in the initialization phase, there are
only three such parameters:
, the drop tolerance used in the incomplete factorization (called dtol in our
code).
maxf il, an integer which controls to overall ll-in (storage) allowed in a given
incomplete factorization.
maxlvl, an integer specifying the maximum number of levels.
(The case corresponds to sparse Gaussian elimina-
tion.) In the solution phase, there are only two additional parameters:
tol, the tolerance used in the convergence test.
maxcg, an integer specifying the maximum number of iterations.
Within our code, all matrices are generally treated within a single, unied frame-
work; e.g., symmetric positive denite, nonsymmetric, and indenite problems generally
do not have specialized options. Besides the control parameters mentioned above,
all information about the matrix is generated from the sparsity pattern and the values
of the nonzeros, as provided in our sparse matrix data structure, a variant of the data
structure introduced in the Yale sparse matrix package [23, 10]. For certain block
matrices, the user may optionally provide a small array containing information about
the block structure.
This input limits the complexity of the code, as well as eliminates parameters
which might be needed to further classify a given matrix. On the other hand, it seems
clear that a specialized solver directed at a specic problem or class of problems, and
making use of this additional knowledge, is likely to outperform our algorithm on
that particular class of problems. Although we do not think our method is provably
\best" for any particular problem, we believe its generality and robustness, coupled
with reasonable computational e-ciency, make it a valuable addition to our collection
of sparse solvers.
The rest of this paper is organized as follows. In section 2, we provide a general
description of our multilevel approach. In section 3, we dene the sparse matrix data
structures used in our code. Our incomplete factorization algorithm is a standard
drop tolerance approach with a few modications for the present application. These
are described in section 4. Our ordering procedure is the minimum degree algorithm.
Once again, our implementation is basically standard with several modications to the
input graph relevant to our application. These are described in section 5. In section
6, we describe the construction of the transfer matrices used in the construction of
the coarse grid correction. Information about the block structure of the matrix, if any
is provided, is used only in the coarsening procedure. This is described in section 7.
Finally, in section 8, we give some numerical illustrations of our method on a variety
of (partial dierential equation) matrices.
2. Matrix formulation. Let A be a large sparse, nonsingular N  N matrix.
We assume that the sparsity pattern of A is symmetric, although the numerical values
need not be. We will begin by describing the basic two-level method for solving
Let B be an N N nonsingular matrix, called the smoother, which gives rise to the
basic iterative method used in the multilevel preconditioner. In our case, B is an
approximate factorization of A, i.e.,
where L is (strict) lower triangular, U is (strict) upper triangular with the same
sparsity pattern as L t , D is diagonal, and P is a permutation matrix.
Given an initial guess x steps of the smoothing procedure produce iterates
The second component of the two-level preconditioner is the coarse grid correction.
Here we assume that the matrix A can be partitioned as
A ff A fc
A cf A cc
where the subscripts f and c denote ne and coarse, respectively. Similar to the
smoother, the partition of A in ne and coarse blocks involves a permutation matrix
P . The ^
coarse grid matrix ^
A is given by
A ff A fc
A cf A cc
I cc
The matrices V cf and W t
N) matrices with identical sparsity patterns;
thus
A has a symmetric sparsity pattern. If A
A.
Let
I cc
In standard multigrid terminology, the matrices ^
W are called restriction and
prolongation, respectively. Given an approximate solution xm to (2.1), the coarse grid
correction produces an iterate xm+1 as follows.
(b Axm );
As is typical of multilevel methods, we dene the two-level preconditioner M
implicitly in terms of the smoother and coarse grid correction. A single cycle takes
an initial guess x 0 to a nal guess x 2m+1 as follows:
Two-Level Preconditioner
are dened using (2.3).
(ii) xm+1 is dened using (2.7).
are dened using (2.3).
The generalization from two-level to multilevel consists of applying recursion to
the solution of the equation ^
r in (2.7). Let ' denote the number of levels in the
recursion.
M(') denote the preconditioner for
A.
Then (2.7) is generalized to
(b Axm );
The general ' level preconditioner M is then dened as follows:
'-Level Preconditioner
directly.
starting from initial guess x 0 , compute x 2m+1 using (iii){(v):
are dened using (2.3).
(iv) xm+1 is dened by (2.8), using iterations of the ' 1 level
scheme for
r to dene ^
M , and with initial guess ^
are dened using (2.3).
The case corresponds to the symmetric V-cycle, while the case
corresponds to the symmetric W-cycle. We note that there are other variants of both
the V-cycle and the W-cycle, as well as other types of multilevel cycling strategies [30].
However, in this work (and in our code) we restrict attention to just the symmetric
V-cycle with presmoothing and postsmoothing iterations.
For the coarse mesh solution our procedure is somewhat nontraditional.
Instead of a direct solution of (2.1), we compute an approximate solution using one
smoothing iteration. We illustrate the practical consequences of this decision in section
8.
If A is symmetric, then so is M , and the '-level preconditioner could be used as a
preconditioner for a symmetric Krylov space method. If A is also positive denite, so
is M , and the standard conjugate gradient method could be used; otherwise the CSCG
method [9], SYMLQ [43], or a similar method could be used. In the nonsymmetric
case, the '-level preconditioner could be used in conjunction with the CSBCG method
[9], GMRES [22], or a similar method.
To complete the denition of the method, we must provide algorithms to
compute the permutation matrix P in (2.2);
compute the incomplete factorization matrix B in (2.2);
compute the ne-coarse partitioning
compute the sparsity patterns and numerical values in the prolongation and
restriction matrices in (2.6).
3. Data structures. Let A be an N N matrix with elements A ij and a symmetric
sparsity structure; that is, both A ij and A ji are treated as nonzero elements
(i.e. stored and processed) if jA diagonal entries A ii are treated as
nonzero regardless of their numerical values.
Our data structure is a modied and generalized version of the data structure
introduced in the (symmetric) Yale sparse matrix package [23]. It is a rowwise version
of the data structure described in [10]. In our scheme, the nonzero entries of A are
stored in a linear array a and accessed through an integer array ja. Let  i be the
number of nonzeros in the strict upper triangular part of row i and set
The array ja is of length N+1+, and the array a is of length N+1+ if A A. If
A t 6= A, then the array a is of length 2. The entries of ja(i), 1  i  N
are pointers dened as follows:
The locations ja(i) to ja(i contain the  i column indices corresponding to
the row i in the strictly upper triangular matrix.
In a similar manner, the array a is dened as follows:
If A t 6= A, then
In words, the diagonal is stored rst, followed by the strict upper triangle stored row-
wise. If A t 6= A, then this is followed by the strict lower triangle stored columnwise.
Since A is structurally symmetric, the column indexes for the upper triangle are identical
to the row indexes for the lower triangle, and hence they need not be duplicated
in storage.
As an example, let
A 11 A 12 A 13 0 0
A 21 A 22 0 A 24 0
A
43 A 44 0
Then
a A11 A22 A33 A44 A55 A12 A13 A24 A34 A35 A21 A31 A42 A43 A53
Diagonal Upper triangle Lower triangle
Although the YSMP data structure was originally devised for sparse direct methods
based on Gaussian elimination, it is also quite natural for iterative methods based
on incomplete triangular decomposition. Because we assume that A has a symmetric
sparsity structure, for many matrix calculations a single indirect address computation
in ja can be used to process both a lower and a upper triangular element in A. For
example, the following procedure computes
procedure mult(N, ja, a, x, y)
end for
for k ja(i) to ja(i
end for
end for
For symmetric matrices, set lmtx 0; umtx 0. Also, may be readily
computed by setting lmtx 0; umtx ja(N
The data structure for storing quite analogous to that
for A. It consists of two arrays, ju and u, corresponding to ja and a, respectively.
The rst entries of ju are pointers as in ja, while entries ju(i) to ju(i
contain column indices of the nonzeros of row i in of U . In the u array, the diagonal
entries of D are stored in the rst N entries. Entry arbitrary. Next, the
nonzero entries of U are stored in correspondence to the column indices in ju. If
the nonzero entries of L follow, stored columnwise.
The data structure we use for the N  ^
W and the ^
are
similar. It consists of an integer array jv and a real array v. The nonzero entries of
are stored rowwise, including the rows of the block I cc . As usual, the rst
entries of jv are pointers; entries jv(i) to jv(i contain column indices for row
W . In the v array, the nonzero entries of ^
are stored rowwise in correspondence
with jv but shifted by N since there is no diagonal part. If ^
W , this is
followed by the nonzeros of ^
stored columnwise.
4. ILU factorization. Our incomplete (L+D)D 1 (D+U) factorization is similar
to the row elimination scheme developed for the symmetric YSMP codes [23, 26].
For simplicity, we begin by discussing a complete factorization and then describe the
modications necessary for the incomplete factorization. Without loss of generality,
assume that the permutation matrix
After k steps of elimination, we have the block factorization
A 11 A 12
A 21 A 22
I
where A 11 is k  k and A 22 is N k  N k. We assume that at this stage, all
the blocks on the right-hand side of (4.1) have been computed except for the Schur
complement S, given by
Our goal for step k + 1 is to compute the rst row and column of S, given by
Because A and (L +D)D 1 (D + U) have symmetric sparsity patterns, and our data
structures take advantage of this symmetry, it is clear that the algorithms for computing
are the same and in practice dier only in the assignments of
shifts for the u and a arrays, analogous to lmtx and umtx in procedure mult. Thus
we will focus on the computation of just . At this point, we also assume that the
array ju has been computed in a so-called symbolic factorization step.
The major substeps are as follows:
1. Copy the rst column of A 22 (stored in the data structures ja and a) into an
expanded work vector z of size N .
2. Find the multipliers given by nonzeros of D 1
3. For each multiplier
using column k of L 21 (i.e.,
4. Copy the nonzeros in z into the data structures ju and u.
In step 1, we need to know the nonzeros of the rst column of A 22 , which is
precisely the information easily accessible in the ja and a data structures. In step
3, we need to know the nonzeros in columns of L 21 , which again is precisely the
information easily available in our data structure. In step 4, we copy a column of
information into the lower triangular portion of the ju and u data structures. Indeed,
the only di-cult aspect of the algorithm is step 2, in which we need to know the
sparsity structure of the rst column of U 12 , information that is not readily available
in the data structure. This is handled in a standard fashion using a dynamic linked
list structure and will not be discussed in detail here.
To generalize this to the incomplete factorization case, we rst observe that the ju
array can be computed concurrently with the numeric factorization simply by creating
a list of the entries of the expanded array z that are updated in step 3. Next, we note
that one may choose which nonzero entries from z to include in the factorization by
choosing which entries to copy to the ju and u data structures in step 4. We do this
through a standard approach using a drop tolerance . In particular, we neglect a
pair of o-diagonal elements if
j. Note D ii has not yet been computed. It is well known that
the ll-in generated through the application of a criterion such as (4.4) is a highly
nonlinear and matrix dependent function of . This is especially problematic in the
present context, since control of the ll-in is necessary in order to control the work
per iteration in the multilevel iteration.
Several authors have explored possibilities of controlling the maximum number
of ll-in elements allowed in each row of the incomplete decomposition [35, 47, 31].
However, for many cases of interest, and in particular for matrices arising from discretizations
of partial dierential equations ordered by the minimum degree algorithm,
most of the ll-in in a complete factorization occurs in the later stages, even if all the
rows initially have about the same number of nonzeros. Thus while it seems advisable
to try to control the total ll-in, one should adaptively decide how to allocate the
ll-in among the rows of the matrix. In our algorithm, in addition to the drop tolerance
, the user provides a parameter maxf il, which species that the total number
of nonzeros in U is not larger than maxf il  N .
Our overall strategy is to compute the incomplete decomposition using the given
drop tolerance. If it fails to meet the given storage bound, we increase the drop
tolerance and begin a new incomplete factorization. We continue in this fashion until
we complete a factorization within the given storage bound. Of course, such repeated
factorizations are computationally expensive, so we developed some heuristics which
allow us to predict a drop tolerance which will satisfy the storage bound.
As the factorization is computed, we make a histogram of the approximate sizes
of all elements that exceed the drop tolerance and are accepted for the factorization.
Let m denote the number of bins in the histogram; our code. Then for
each pair of accepted o-diagonal elements, we nd the largest k 2 [1; m] such that
Here  > 1 our code). The histogram is realized as an integer array h of
size m, where h ' is the number of accepted elements that exceeded the drop tolerance
by factors between  ' 1 and  ' for 1  '  m 1; hm contains the number of
accepted elements exceeding the drop tolerance by  m 1 . If the factorization reaches
the storage bound, we continue the factorization but allow no further ll-in. However,
we continue to compute the histogram based on (4.5), proling the elements we would
have accepted had space been available. Then using the histogram, we predict a
new value of  such that the total number of elements accepted for U is no larger
than maxf il  N=. Such a prediction of course cannot be guaranteed, since the
sizes and numbers of ll-in elements depend in a complicated fashion on the specic
history of the incomplete factorization process; indeed, the histogram cannot even
completely prole the remainder of the factorization with the existing drop tolerance,
since elements that would have been accepted could introduce additional ll-in at
later stages of the calculation as well as in
uence the sizes of elements computed at
later stages of the factorization. In our implementation, the factor  varies between
depending on how severely the storage bound was exceeded. Its
purpose is to introduce some conservative bias into the prediction with the goal that
the actual ll-in accepted should not exceed maxf il  N .
Finally, we note that there is no comprehensive theory regarding the stability
of incomplete triangular decompositions. For certain classes of matrices (e.g., M-matrices
and H-matrices), the existence of certain incomplete factorizations has been
proved [39, 25, 24, 40, 51]. However, in the general case, with potentially indenite
and/or highly nonsymmetric matrices, one must contend in a practical way with
the possibility of failure or near failure of the factorization. A common approach
is to add a diagonal matrix, often a multiple of the identity, to A and compute an
incomplete factorization of the shifted matrix. One might also try to incorporate
some form of diagonal pivoting; partial or complete pivoting could potentially destroy
the symmetric sparsity pattern of the matrix. However, any sort of pivoting greatly
increases the complexity of the implementation, since the simple but essentially static
data structures ja, a, ju, and u are not appropriate for such an environment.
Our philosophy here is to simply accept occasional failures and continue with
the factorization. Our ordering procedure contains some heuristics directed towards
avoiding or at least minimizing the possibility of failures. And when they do occur,
failures often corrupt only a low dimensional subspace, so a Krylov space method
such as conjugate gradients can compensate for such corruption with only a few extra
iterations. In our implementation, a failure is revealed by some diagonal entries in D
becoming close to zero. O-diagonal elements L ji and U ij are multiplied by D 1
ii , and
the solution of (L multiplication by D 1
ii . For
purposes of calculating the factorization and solution, the value of D 1
ii is modied
near zero as follows:
1=D ii for jD ii j > ,
Here  is a small constant; in our implementation,  is the machine
epsilon. Although many failures could render the preconditioner well-dened but
essentially useless, in practice we have noted that D 1
ii is rarely modied for the large
class of nite element matrices which are the main target of our procedure.
5. Ordering. To compute the permutation matrix P in (2.2), we use the well-known
minimum degree algorithm [45, 26]. Intuitively, if one is computing an incomplete
factorization, an ordering which tends to minimize the ll-in in a complete
factorization should tend to minimize the error
For particular classes of matrices, specialized ordering schemes have been developed
[34, 15, 37, 36]. For example, for matrices arising from convection dominated prob-
lems, ordering along the
ow direction has been used with great success. However, in
this general setting, we prefer to use just one strategy for all matrices. This reduces
the complexity of the implementation and avoids the problem of developing heuristics
to decide among various ordering possibilities. We remark that for convection
dominated problems, minimum degree orderings perform comparably well to the specialized
ones, provided some (modest) ll-in is allowed in the incomplete factorization.
For us, this seems to be a reasonable compromise.
Our minimum degree ordering is a standard implementation, using the quotient
graph model [26] and other standard enhancements. A description of the graph of
the matrix is the main required input. Without going into detail, this is essentially
a small variant of the basic ja data structure used to store the matrix A. We will
denote this modied data structure as jc. Instead of storing only column indices for
the strict upper triangle as in ja, entries jc(i) to jc(i of the jc data structure
contain column indices for all o-diagonal entries for row i of the matrix A.
We have implemented two small enhancements to the minimum degree ordering;
as a practical matter, both involve changes to the input graph data structure jc
that is provided to the minimum degree code. First, we have implemented a drop
tolerance similar to that used in the the factorization. In particular the edge in the
graph corresponding to o-diagonal entries A ij and A ji is not included in the jc data
structure if
jA jj A ii j: (5.1)
This excludes many entries which are likely to be dropped in the subsequent incomplete
factorization and hopefully will result in an ordering that tends to minimize the
ll-in created by the edges that are kept.
The second modication involves some modest a priori diagonal pivoting designed
to minimize the number failures (near zero diagonal elements) in the subsequent
factorization. We rst remark that pivoting or other procedures based on the values
of the matrix elements (which can be viewed as weights on graph edges and nodes)
would destroy many of the enhancements which allow the minimum degree algorithm
to run in almost linear time. Our modication is best explained in the context of a
simple 2  2 example. Let
b a
with a; b; c 6= 0. Clearly, A is nonsingular, but the complete triangular factorization
of A does not exist. However,
a b
a 0
c bc=a
a b
Now suppose that A ii  0, A jj these four elements form a
submatrix of the form described above, and it seems an incomplete factorization of
A is less likely to fail if the P is chosen such that vertex j is ordered before vertex i.
This is done as follows: for each i such that A ii  0, we determine a corresponding
j such that A jj there is more than one choice, we choose the one for
which jA ij A ji =A jj j is maximized. To ensure that vertex i is ordered after vertex j,
we replace the sparsity pattern for the o-diagonal entries for row (column) i with
the union of those for rows (columns) i and j. If we denote the set of column indices
for row i in the jc array as adj(i), then
Although the sets adj(i) and adj(j) are modied at various stages, it is well known
that (5.3) is maintained throughout the minimum degree ordering process [26], so
that at every step of the ordering process deg(j)  deg(i), where deg(i) is the degree
of vertex i. As long as deg(j) < deg(i), vertex j will be ordered before vertex i by
the minimum degree algorithm. On the other hand, if at some stage
of the ordering process, it remains so thereafter, and (5.3) becomes
In words, i and j become so-called equivalent vertices and will be eliminated at the
same time by the minimum degree algorithm (see [26] for details). Since the minimum
degree algorithm sees these vertices as equivalent, they will be ordered in an arbitrary
fashion when eliminated from the graph. Thus, as a simple postprocessing step, we
must scan the ordering provided by the minimum degree algorithm and exchange the
order of rows i and j if i was ordered rst. Any such exchanges result in a new
minimum degree ordering which is completely equivalent, in terms of ll-in, to the
the original.
For many types of nite element matrices (e.g., the indenite matrices arising
from Helmholtz equations), this a priori scheme is useless because none of the diagonal
entries of A is close to zero. However, this type of problem is likely to produce only
isolated small diagonal entries in the factorization process, if it produces any at all.
On the other hand, other classes of nite element matrices, notably those arising
in from mixed methods, Stokes equations, and other saddle-point-like formulations,
have many diagonal entries that are small or zero. In such cases, the a priori diagonal
pivoting strategy can make a substantial dierence and greatly reduce the numbers
of failures in the incomplete triangular decomposition.
6. Computing the transfer matrices. There are three major tasks in computing
the prolongation and restriction matrices ^
W of (2.6). First, one must
determine the sparsity structure of these matrices; this involves choosing which unknowns
are coarse and which are ne. This reduces to determining the permutation
P of (2.4). Second, one must determine how coarse and ne unknowns are
related, the so-called parent-child relations [49]. This involves computing the sparsity
patterns for the matrices V cf and W fc . Third, one must compute the numerical values
for these matrices, the so-called interpolation coe-cients [50].
There are many existing algorithms for coarsening graphs. For matrices arising
from discretizations of partial dierential equations, often the sparsity of the matrix A
is related in some way to the underlying grid, and the problem of coarsening the graph
of the matrix A can be formulated in terms of coarsening the grid. Some examples are
given in [14, 13, 17, 18, 46, 12, 49]. In this case, one has the geometry of the grid to
serve as an aid in developing and analyzing the coarsening procedure. There are also
more general graph coarsening algorithms [32, 33, 19], often used to partition problems
for parallel computation. Here our coarsening scheme is based upon another well-known
sparse matrix ordering technique, the reverse Cuthill{McKee algorithm. This
ordering tends to yield reordered matrices with minimal bandwidth and is widely used
with generalized band elimination algorithms [26]. We now assume that the graph
has been ordered in this fashion and that a jc data structure representing the graph
in this ordering is available. Our coarsening procedure is just a simple postprocessing
step of the basic ordering routine, in which the N vertices of graph are marked as
COARSE or F INE.
procedure coarsen(N, jc, type)
end for
for j jc(i) to jc(i
end for
end for
This postprocessing step, coupled with the the reverse Cuthill{McKee algorithm,
is quite similar to a greedy algorithm for computing maximal independent sets using
breadth-rst search. Under this procedure, all coarse vertices are surrounded only by
ne vertices. This implies that the matrix A cc in (2.4) is a diagonal matrix. For the
sparsity patterns of matrices arising from discretizations of scalar partial dierential
equations in two space dimensions, the number of coarse unknowns ^
N is typically on
the order of N=4 to N=5. Matrices with more nonzeros per row tend to have smaller
values of ^
N . To dene the parents of a coarse vertex, we take all the connections of
the vertex to other ne vertices; that is, the sparsity structure of V cf in (2.5) is the
same as that of the block A cf .
In our present code, we pick V cf and W fc according to the formulae
ff A fc ;
ff
~
Here D ff is a diagonal matrix with diagonal entries equal to those of A ff . In this
sense, the nonzero entries in V cf and W fc are chosen as multipliers in Gaussian elimi-
nation. The nonnegative diagonal matrices R ff and ~
R ff are chosen such that nonzero
rows of W fc and columns of V cf , respectively, have unit norms in ' 1 .
Finally, the coarsened matrix ^
A of (2.5) is \sparsied" using the drop tolerance
and a criterion like (5.1) to remove small o-diagonal elements. Empirically, applying
a drop tolerance to ^
A at the end of the coarsening procedure has proved more e-cient,
and more eective, than trying to independently sparsify its constituent matrices. If
the number of o-diagonal elements in the upper triangle exceeds maxf il  ^
N , the
drop tolerance is modied in a fashion similar to the incomplete factorization. The
o-diagonal elements are proled by a procedure similar to that for the incomplete fac-
torization, but in this case the resulting histogram is exact. Based on this histogram,
a new drop tolerance is computed, and (5.1) is applied to produce a coarsened matrix
satisfying the storage bound.
7. Block matrices. Our algorithm provides a simple but limited functionality
for handling block matrices. Suppose that the N N matrix A has the K K block
structure
A =B @
where subscripts for A ij are block indices and the diagonal blocks A jj are square
matrices. Suppose A jj is of order
The matrix A is stored in the usual ja and a data structures as described in
section 3 with no reference to the block structure. A small additional integer array ib
of size K + 1 is used to dene the block boundaries as follows:
In words, integers in the range ib(j) to ib(j inclusive, comprise the index set
associated with block A jj . Note that ib(K
This block information plays a role only in the coarsening algorithm. First, the
reverse Cuthill{McKee algorithm described in section 6 is applied to the block diagonal
matrix
A =B @
A 11
AKKC A (7.2)
rather than A. As a practical matter, this involves discarding graph edges connecting
vertices of dierent blocks in the construction of the graph array jc used as input.
Such edges are straightforward to determine from the information provided in the ib
array. The coarsening algorithm applied to the graph of
A produces output equivalent
to the application of the procedure independently to each diagonal block of
A. As
a consequence, the restriction and prolongation matrices automatically inherit the
block structure A. In particular,
VKKC A and ^
W jj are are rectangular matrices
having the structure of (2.6) that would have resulted from the application of the
algorithm independently to A jj . However, like the matrix A,
are stored
in the standard jv and v data structures described in section 3 without reference to
their block structures.
The complete matrix A is used in the construction of the coarsened matrix ^
A of
(2.5). However, because of (7.1) and (7.3)
A 1K
A also automatically inherits the K K block structure of A. It is not necessary
for the procedure forming ^
A to have any knowledge of its block structure, as this
block structure can be computed a priori by the graph coarsening procedure. Like
A is stored in standard ja and a data structures without reference to its block
structure. Since the blocks of A have arbitrary order, and are essentially coarsened
independently, it is likely that eventually some of the ^
That is, certain blocks
may cease to exist on coarse levels. Since the block information is used only to discard
certain edges in the construction of the graph array jc, \00" diagonal blocks present
no di-culty.
8. Numerical experiments. In this section, we present a few numerical illus-
trations. In our rst sequence of experiments, we consider several matrices loosely
based on the classical case of 5-point centered nite dierence approximations to u
on a uniform square mesh. Dirichlet boundary conditions are imposed. This leads to
the n  n block tridiagonal system
I
I T I
I T I
I
with T the n  n tridiagonal matrix
This is a simple test problem easily solved by standard multigrid methods. In contrast
to this example we also consider the block tridiagonal system
Both A and
A have the same eigenvectors and the same eigenvalues, although the
association of eigenvectors and eigenvalues are reversed in the case of
A. That is,
the so-called smooth eigenvectors are associated with large eigenvalues, while rough
eigenvectors are associated with smaller eigenvalues. Although
A does not arise naturally
in the context of numerical discretizations of partial dierential equations, it is
of interest because it dees much of the conventional wisdom for multigrid methods.
Third, we consider block 3  3 systems of the form
where A is the discrete Laplacian and D is a symmetric positive denite \stabilization"
matrix with a sparsity pattern similar to A. However, the nonzeros in D are of size
compared to size O(1) nonzero elements in A. C x and C y also have sparsity
patterns similar to that of A, but these matrices are nonsymmetric and their nonzero
entries are of size O(h). Such matrices arise in stabilized discretizations of the Stokes
equations. One third of the eigenvalues of S are negative, so S is quite indenite.
In addition to the ja and a arrays, for the matrix S we also provided an ib array as
described in section 7 to dene its 33 block structure. We emphasize again that this
block information is used only in the computation of the graph input to the coarsening
procedure and is not involved in any aspect of the incomplete factorization smoothing
procedure. With many small diagonal elements, this class of matrices provides a good
test of the a priori pivoting strategy used in conjunction with the minimum degree
ordering.
In

Table

8.1, Levels refers to the number of levels used in the calculation. In our
implementation the parameter maxlvl, which limits the number of levels allowed, was
set su-ciently large that it had no eect on the computation. The drop tolerance
was set to matrices. The ll-in control parameter maxf il was set
su-ciently large that it had no eect on the computation. The initial guess for all
problems was x
In

Table

8.1, the parameter Digits refers to
In these experiments, we asked for six digits of accuracy. The column labeled Cycles
indicates the number of multigrid cycles (accelerated by CSCG) that were used to
achieve the indicated number of digits. Finally, the last two columns, labeled Init. and
Solve, record the CPU time, measured in seconds, for the initialization and solution
phases of the algorithm, respectively. Initialization includes all the orderings, incomplete
factorizations, and computation of transfer matrices used in the multigraph
preconditioner. Solution includes the time to solve (2.1) to at least six digits given
the preconditioner. These experiments were run on an SGI Octane R10000 250mhz,
using double precision arithmetic and the f90 compiler.
In analyzing these results, it is clear that our procedure does reasonably well
on all three classes of matrices. Although it appears that the rate of convergence
is not independent of N , it seems apparent that the work is growing no faster than
logarithmically. CPU times for larger vales of N are aected by cache performance
as well as the slightly larger number of cycles.
For the highly indenite Stokes matrices S, it is important to also note the ro-
bustness, that the procedure solved all of the problems. With more nonzeros per row
on average, the incomplete factorization was more expensive to compute than for the
other cases. This is re
ected in relatively larger initialization and solve times.
In our next experiment, we illustrate the eect of the parameters maxlvl and .
For the matrix A with we solved the problem for
and 1  maxlvl  7. We terminated the iteration when the solution had six digits,
Performance comparison.
Digits Cycles Init. Solve
Discrete Laplacian A,
Stokes matrix S,
as measured by (8.1). We also provide the total storage for the ja and ju arrays
for all matrices, measured in thousands of entries. Since the matrices are symmetric,
this is also the total (
oating point) storage for all matrices A and approximate LDU
factorizations.
Here we see that our method behaves in a very predictable way. In particular,
decreasing the drop tolerance or increasing the number of levels improves the convergence
behavior of the method. On the other hand, the timings do not always follow
the same trend. For example, for the case increasing the number of levels
from decreases the number of cycles but increases the time.
This is because for our method defaults to the standard conjugate gradient
iteration with the incomplete factorization preconditioner. When maxlvl > 1,
one presmoothing and one postsmoothing step are used for the largest matrix. With
the additional cost of the recursion, the overall cost of the preconditioner is more than
double the cost for the case
We also note that, unlike the classical multigrid method, where the coarsest matrix
is solved exactly, in our code we have chosen to approximately solve the coarsest
system using just one smoothing iteration using the incomplete factorization. When
the maximum number of levels are used, as in Table 8.1, the smallest system is
typically 1  1 or 2  2, and this is an irrelevant remark. However, in the case
of

Table

8.2, the fact that the smallest system is not solved exactly signicantly
in
uences the overall rate of convergence. This is why, unlike methods where the
coarsest system is solved exactly, increasing the number of levels tends to improve
the rate of convergence. In the case the coarsest matrix had an exact LDU
factorization for the case (because the matrix itself was nearly diagonal),
and setting maxlvl > 5 did not increase the number of levels. The cases
Dependence of convergence of  and maxlvl, discrete Laplacian A,
maxlvl Digits Cycles Init. Solve
3 6.1 96 13.2 116.9 1077 1119
6 { { { { { {
7 { { { { { {
6.1 56 12.1 64.9 878 2106
6.1 22 16.6 31.7 878 3649
and used a maximum of 10 and 9 levels, respectively, but the results did not
change signicantly from the case 7.
We also include in Table 8.2 the case
ination. (In fact, our code uses jjAjj as the drop tolerance when the user species
to avoid dividing by zero.) Here we see that Gaussian elimination is reasonably
competitive on this problem. However, we generally expect the initialization cost for
to grow like O(N 3=2 ). For we expect the solution times to
grow like O(N p ), p > 1. For the best multilevel choices, we expect both initialization
and solution times to behave like O(N) O(N log N ).
In our nal series of tests, we study the convergence of the method for a suite of
test problems generated from the nite element code PLTMG [8]. These example
problems were presented in our earlier work [11], where a more complete description
of the problems, as well as numerical results for our hierarchical basis multigraph
method and the classical AMG algorithm of Ruge and Stuben [46], can be found.
As a group, the problems feature highly nonuniform, adaptively generated meshes,
relatively complicated geometry, and a variety of dierential operators. For each test
case, both the sparse matrix and the right-hand side were saved in a le to serve as
input for the iterative solvers. A short description of each test problem is given below.
Problem Superior. This problem is a simple Poisson equation
with homogeneous Dirichlet boundary conditions on a domain in the shape of Lake
Superior. This is the classical problem on a fairly complicated domain. The solution
is generally very smooth but has some boundary singularities.
Problem Hole. This problem features discontinuous, anisotropic coe-cients. The
overall domain is the region between two concentric circles, but this domain is divided
into three subregions. On the inner region, the problem is
with In the middle region, the equation is
and in the outer region the equation is
Homogeneous Dirichlet boundary conditions are imposed on the inner (hole) bound-
ary, homogeneous Neumann conditions on the outer boundary, and the natural continuity
conditions on the internal interfaces. While the solution is also relatively
smooth, singularities exist at the internal interfaces.
Problem Texas. This is an indenite Helmholtz equation
posed in a region shaped like the state of Texas. Homogeneous Dirichlet boundary
conditions are imposed. The length scales of this domain are roughly 16  16, so this
problem is fairly indenite.
Problem UCSD. This is a simple constant coe-cient convection-diusion equation
r  (ru
posed on a domain in the shape of the UCSD logo. Homogeneous
Dirichlet boundary conditions are imposed. Boundary layers are formed at the bottom
of the region and the top of various obstacles.
Problems Jcn 0 and Jcn 180. The next two problems are solutions of the current
continuity equation taken from semiconductor device modeling. This equation is a
convection-diusion equation of the form
r  (ru
in most of the rectangular domain. However, in a curved band in the interior of
the domain, jj  10 4 and is directed radially. Dirichlet boundary conditions
and are imposed along the bottom boundary and along a short segment on
the upper left boundary, respectively. Homogeneous Neumann boundary conditions
are specied elsewhere. The solutions vary exponentially across the domain which is
typical of semiconductor problems.
In the rst problem, Jcn 0, the convective term is chosen so the device is forward
biased. In this case, a sharp internal layer develops along the top interface boundary.
In the second problem, Jcn 180, the sign of the convective term is reversed, resulting
in two sharp internal layers along both interface boundaries.
We summarize the results in Table 8.3. As before, perhaps the most important
point is that the method solved all of the problems. While convergence rates are not
independent of h, once again the growth appears to be at worst logarithmic.
Below we make some additional remarks.

Table
Performance comparison.
N Levels Digits Cycles Init. Solve
20k 9 7.3 5 1.4e 0 9.4e-1
Hole,
Jcn 0,
Jcn
For all problems, decreasing the drop tolerance will tend to increase the effectiveness
of the preconditioner, although it generally will also make the
preconditioner more costly to apply. Thus one might optimize the selection
of the drop tolerance to minimize the decreasing number of cycles against the
increasing cost per cycle. In these experiments, we did not try such systematic
optimization, but we did adjust the drop tolerance in a crude way such
that more di-cult problems performed in a fashion similar to the easy ones.
Problem Texas is by far the most di-cult in this test suite. While we set
the problem with order 80k was the only one which came close
to achieving this storage limit. Most were well below this limit, and many
averaged less than 10 nonzeros per row in L and U factors.
For the nonsymmetric problems the CSBCG method is used for acceleration.
Since the CSBCG requires the solution of a conjugate system with A t , two
matrix multiplies and two preconditioning steps are required for each itera-
tion. As noted in section 3, with our data structures, applying a transposed
matrix and preconditioner costs the same as applying the original matrix or
preconditioner. Since these are the dominant costs in the CSBCG methods,
the cost per cycle is approximately double that for an equivalent symmetric
system.



--R

On eigenvalue estimates for block incomplete factorization methods

The algebraic multilevel iteration methods - theory and applications
Stabilization of algebraic multilevel iteration methods
Algebraic multilevel preconditioning methods I

A class of hybrid algebraic multilevel preconditioning methods
PLTMG: A Software Package for Solving Elliptic Partial Di
An analysis of the composite step biconjugate gradient method
General sparse elimination requires no permanent integer storage
The incomplete factorization multigraph algorithm

The hierarchical basis multigrid method and incomplete LU decomposi- tion

Orderings for incomplete factorization preconditioning of nonsymmetric problems
Towards algebraic multigrid for elliptic problems of second order


Boundary treatments for multilevel methods on unstructured meshes

Black box multigrid
Variational iterative methods for non-symmetric systems of linear equations
Algorithms and data structures for sparse symmetric Gaussian elimination
A stability analysis of incomplete LU factorizations
Algebraic analysis of the hierarchical basis preconditioner
Computer Solution of Large Sparse Positive De
Incomplete block factorization preconditioning for linear systems arising in the numerical solution of the helmholtz equation

An algebraic hierarchical basis preconditioner

Incomplete Decompositions - Theory
Analysis of multilevel graph partitioning

Ordering techniques for convection dominated problems on unstructured three dimensional grids

Ordering strategies for modi

Energy optimization of algebraic multigrid bases
An analysis of the robustness of some incomplete factorizations
On the stability of the incomplete LU-factorizations and characterizations of H-matrices
Using approximate inverses in algebraic multilevel methods

Solution of sparse inde
A multigrid method based on incomplete Gaussian elimination
A graph theoretic study of the numeric solution of sparse positive de

ILUT: a dual threshold incomplete LU factorization
Convergence of algebraic multigrid based on smoothed aggregation
Introduction to algebraic multigrid
An energy-minimizing interpolation for robust multi-grid methods
On the robustness of ILU smoothing
--TR

--CTR
Randolph E. Bank, Compatible coarsening in the multigraph algorithm, Advances in Engineering Software, v.38 n.5, p.287-294, May, 2007
Gh. Juncu , E. Mosekilde , C. Popa, Numerical experiments with MG continuation algorithms, Applied Numerical Mathematics, v.56 n.6, p.844-861, June 2006
J. S. Ovall, Hierarchical matrix techniques for a domain decomposition algorithm, Computing, v.80 n.4, p.287-297, September 2007
Michele Benzi, Preconditioning techniques for large linear systems: a survey, Journal of Computational Physics, v.182 n.2, p.418-477, November 2002
algebraic multigrid;multigraph methods;incomplete LU factorization