The Incomplete Factorization Multigraph Algorithm.
We present a new family of multigraph algorithms, ILU-MG, based upon an incomplete sparse matrix factorization using a particular ordering and allowing a limited amount of fill-in. While much of the motivation for multigraph comes from multigrid ideas, ILU-MG is distinctly different from algebraic multilevel methods. The graph of the sparse matrix A is recursively coarsened by eliminating vertices using a graph model similar to Gaussian elimination. Incomplete factorizations are obtained  by allowing only the fill-in generated by the vertex parents associated with each vertex. Multigraph is numerically compared with algebraic multigrid on some examples arising from discretizations of partial differential equations on unstructured grids.
Introduction
. In this paper, we present a new family of multigraph algo-
rithms, ILU-MG, based upon an incomplete sparse matrix factorization using a carefully
designed ordering and allowing a limited amount of fill-in. While in this paper we
focus primarily on systems of linear equations arising from discretizations of partial
differential equations, the method can be formally applied to general sparse matri-
ces. For any particular problem or class of problems, it seems likely that specialized
methods making use of the particular features of that problem will outperform any
multigraph algorithm. However, the goal of the ILU-MG algorithm is to provide a
general and robust iterative solver for many different systems of linear equations.
While this goal may not yet be achieved in this first version, our hope and expectation
is that the multigraph algorithm will eventually provide reasonably good rates
of convergence for many classes of problems, while requiring only minimal input.
Algebraic approaches to multilevel methods have enjoyed a long history, beginning
with the algebraic multigrid (AMG) methods of Brandt, McCormick, and Ruge [13,
14] Ruge and St-uben [26] and the black box multigrid method of Dendy [15]. More
recent work can be found in [1, 3, 4, 12, 20, 19, 17], as well as many contributions in
[2]. Our work grew out of the grid coarsening schemes developed in [10, 11] and the
corresponding hierarchical basis iterations, HBMG. While much of our motivation
comes from these multigrid ideas, ILU-MG is fundamentally an incomplete sparse
matrix factorization.
The multigraph method resembles the approach of classical sparse Gaussian elimi-
nation. The graph of the stiffness matrix A is recursively coarsened by first eliminating
a node with all its adjacent edges and then adding only a partial set of fill-in edges
corresponding to the vertex parents of that node. Thus, there is no concept of levels,
and more important, no coarse grid on which the problem must be solved exactly.
However, the generality of this approach leaves open the possibility of introducing
both levels and coarse graphs. We plan to study these alternatives as possible ways
to improve the convergence behavior of the basic method. See [9] for some preliminary
results in this direction. The size of each subsequent graph is controlled by
monitoring the amount of numerical fill-in produced. Since each node is eliminated
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.
based upon graph considerations, no attempt is made to preserve the integrity of the
grid. Indeed, this information is not even provided.
The rest of this paper is organized as follows. In section 2, we provide a graph
theoretic interpretation of the construction of hierarchical bases and its relation to
sparse incomplete factorizations. The connection between multigrid and multigraph
is described in terms of linear algebra in section 3. In section 4, the implementation
of our method is discussed in some detail. In particular, the ordering strategy and
the incomplete factorization procedure are described. Finally, in section 5, we compare
multigraph with AMG on some examples arising from discretizations of partial
differential equations on unstructured grids.
2. Graph theoretical aspects. In this section we discuss the relation between
the construction of a hierarchical basis and sparse incomplete LU (ILU) factorization,
within the context of graph theory. We first consider standard Gaussian elimination
and classical ILU factorization from a graph theoretical point of view, and then
develop a graph elimination model for hierarchical basis methods on sequences of
nested meshes. These models can be interpreted as special ILU decompositions which
generalize to the case of general graphs.
We begin with a few standard definitions; the interested reader is referred to
Rose [25] or George and Liu [18] for a more complete introduction. Corresponding to
a sparse n \Theta n matrix A with symmetric sparsity pattern (i.e., A ij 6= 0 if and only if
A ji 6= 0), let G(V; E) be the graph that consists of a set of n ordered vertices
set of edges E such that the edge (connecting vertices v i and
j. The edges in the graph G correspond to the nonzero
off-diagonal entries of A. If A is the stiffness matrix for the space of continuous
piecewise linear polynomials represented in the standard nodal basis, the graph G
is just the underlying triangulation of the domain (with minor modifications due to
Dirichlet boundary conditions). For vertex v i the set of adjacent vertices adj(v i ) is
defined by
The degree of a vertex deg(v i ) is just the size of the set adj(v i ). A clique C ' V is a
set of vertices which are all pairwise connected; that is, v
With a proper ordering of the vertices, a clique corresponds to a dense submatrix of
A. In graph theoretic terms, a single step of Gaussian elimination transforms G(V; E)
to a new graph G 0
1. Eliminate vertex v i and all its incident edges from G. Set g.
Denote the resulting set of edges
2. Create a set F of fill-in edges as follows: for each distinct pair v
in G, add the edge e jk to F if not already present in E 1 . Set
Since the values of matrix entries are not involved, this model cannot take into
account the occurrence of so-called accidental zeros. The graph elimination process is
illustrated in Figure 1. Note that the set adj(v) in G becomes a clique in G 0 . Within
this framework, the classical ILU factorization is one in which no fill-in edges are
allowed, i.e., F j ;. This forces the matrix A 0 corresponding to the new graph G 0 to
have the same sparsity structure as the corresponding submatrix of A. The graph G 0
would then correspond to the center picture in Figure 1.
The concept of vertex parents is first introduced to allow HBMG to be interpreted
as a generalized ILU procedure. We will begin with the case of two nested meshes
where the fine mesh is a uniform refinement of a coarse mesh, generated by pairwise
\Phi \Phi \Phi \Phi \Phi
\Phi\Omega \Omega \Omega \Omega
\Omega \Omega \Omega \Omega \Phi \Phi \Phi
Fig. 1. A fine grid vertex v is eliminated by classical Gaussian elimination. The original mesh
is shown on the left. We remove v and its incident edges (middle), and add fill-in edges (right).
connecting the midpoints of the coarse grid edges in the usual way [7, 27, 21]. Here
we can make the direct sum decomposition c is the set of coarse
grid vertices and V f is the set of fine grid vertices (those not in V c ). For each vertex
there is a unique pair of vertex parents v c such that v i is the midpoint
of the edge connecting v j and v k (v
We now view HBMG as an ILU algorithm in which only selected fill-in edges are
allowed. In this algorithm, the vertices in the set V f are sequentially eliminated as
follows:
1. Eliminate vertex its incident edges from G. Set g.
Denote the resulting set of edges
2. Let v the parents of v i . Create a set F of fill-in edges
of the form e already
present in E 1 . Set
In other words, the classical HBMG algorithm adds the subset of fill-in edges from
Gaussian elimination in which one of the vertices is a vertex parent.
An even more simple possibility is to use just one vertex parent as an elimination
strategy. Although this does not correspond to the classical case of HBMG, it has
been studied in a different context as a partitioning scheme for general graphs [23, 24].
The elimination algorithm is similar to the case of two parents:
1. Eliminate vertex its incident edges from G. Set g.
Denote the resulting set of edges
2. Let v j 2 adj(v i ) denote the parent of v i . Create a set F of fill-in edges of
the form e already present in E 1 . Set
However, in this case generally fewer fill-in edges are added. When the initial graph
G is a finite element triangulation (or a tetrahedral mesh in three space dimensions),
the graph G 0 remains a finite element triangulation. This property can be maintained
at all steps of the elimination process through a careful selection of parents. Both one
and two vertex parent eliminations are illustrated in Figure 2.
Another straightforward extension allows for the possibility of more than two
vertex parents, either a fixed or variable number. The limiting case, allowing all
vertices in adj(v) to be parents of v, will result in classical Gaussian elimination.
These extensions can be interpreted, in the multigrid framework, as various multipoint
interpolation schemes. In this work, we consider only the one- and two-vertex parent
cases.
Let the triangulation T f be the graph for the original stiffness matrix A represented
in the standard nodal basis. For either one or two parents, after all the
\Phi\Omega \Omega \Omega \Omega
\Omega \Omega \Omega \Omega \Phi \Phi \Phi \Phi \Phi
\Phi\Omega \Omega \Omega \Omega
Fig. 2. A fine grid vertex v is eliminated by HBMG ILU elimination. The original mesh is
shown on the left. The fill-in pattern using two parents (v 1 and v 2 ) is shown in the middle, while
the fill-in pattern using the single parent v 1
is shown on the right.
vertices in V f are eliminated, the resulting graph is just the coarse grid triangulation
T c . However, the numerical values of the matrix elements are generally different for
the two cases. For the special case of a sequence of uniformly refined meshes, the
total number of edges in the filled in graph can be estimated. This will serve as a
guide to the amount of memory necessary to store the incomplete LU factorization
using typical sparse matrix storage schemes, e.g., [16, 18, 8]. The elimination process
is illustrated for the case of two parents in Figure 3.
A
A
A
A
A
A
A
A
A
A
A
A
A
A
A
A
A
A
A
A
A
A
A
Fig. 3. HBMG=ILU in the classical case, using two parents. The fine grid vertices are
eliminated in the order v 1 , v 2 , v 3 . A total of two fill-in edges are added in the interior of the coarse
grid triangle during the elimination.
The original graph has approximately 3n edges and 2n triangles. Each quartet
of four triangles generates two fill-in edges which are not part of the coarse grid
triangulation. Thus the total number of edges not in the coarse grid is approximately
4n. If we repeat recursively for the coarser grids, the geometric sum
4n
ae
oe
16nis generated. Thus there are approximately 16n=3 edges in the filled in graph. In
the case of classical HBMG, the method is implemented as a block iteration, with the
fill-in in the off-diagonal blocks not explicitly stored. The number of edges actually
stored in the sparse data structure is approximately shown in [7]. In
our present study, we consider only point iterations with explicit storage of all edges.
The elimination process for the case of one parent is shown in Figure 4.
The main difference in the case of one parent is that now only one fill-in edge is
associated with each quartet of triangles on the fine grid. Thus we have a total of
approximately edges not in the coarse grid and expect approximately
7nae
oe
A
A
A
A
A
A
A
A
A
A
A
A
A
A
A
A
A
A
A
A
A
A
A
A
A
A
A
A
A
A
A
A
A
Fig. 4. HBMG ILU using one parent. The fine grid vertices are eliminated in the order v 1
v 3 . There are eight generic scenarios for the process, but all result in one fill-in edge in the interior
of the coarse grid triangle.
total edges, somewhat less than the case of two parents. This is also less than the
bound of 6n obtained for general matching strategies on two-dimensional triangulations
[23]. In fact, matching strategies when applied to general graphs produce fill-in
bounded by O(jE is the number of edges in the original graph. We
know of no similar bound for the case of two parents, and have empirically observed
very rapid growth in fill-in for some matching strategies which made little effort to
control the number of fill-in edges.
The multigraph algorithm generalizes the concept of vertex parents to arbitrary
graphs. The main problem is to determine reasonable vertex parents for each vertex
to be eliminated. Once this is done, the elimination/unrefinement/coarsening is
performed on the graph exactly as in the case of nested meshes.
3. The multigraph algorithm. The salient features of the multigraph algorithm
and its connection to HBMG can be illustrated by considering the following
simple example. Let A be an n \Theta n sparse matrix arising from the discretization of a
partial differential equation, assembled using the standard nodal basis. We partition
A as
(1)
where d 6= 0 is the diagonal matrix element for vertex v 1 , r and c are vectors of order
is an (n \Gamma 1) \Theta (n \Gamma 1) matrix. The hierarchical basis multigrid method
is based on a change of basis, from the nodal basis to a hierarchical basis. In the
present context, the first step of this process involves forming the matrix
I
I
(2)
The vectors ' and u are sparse with nonzeros determined by the vertex parents of
vertex
For classical HBMG, v 1 is a vertex associated with the refinement of an edge with
endpoints g. In this case,
is the usual unit vector (jth column of the identity matrix I Here the
restriction and interpolation operators are chosen geometrically, reflecting the fact
that
One possible generalization is to choose the vectors ' and u based on an incomplete
factorization of the matrix A. In particular, let v i and v j be the parents of vertex v 1
and
The scalars ' are simply the multipliers of an incomplete LU factoriza-
tion. Note that the above equations are well defined, independent of the geometry of
the mesh, and they do not require that v 1 results from the refinement of an edge e ij .
In a similar manner, for the case of one vertex parent we have
The sparsity patterns of the vectors c and c + d' are the same, since the nonzeros
in ' are a subset of the nonzeros in c. 1 However, the sparsity patterns of B and
generally are not the same. The matrix 'du t typically
creates some fill-in in the rows and columns corresponding to the vertex parents.
These are precisely the fill-in edges illustrated in Figure 2. Also note that
From this identity, we can see that this elimination step can also be viewed as forming
a rank one perturbation of the exact Schur complement
In its standard formulation, the next step of this hierarchical decomposition is
a transformation of the same form applied to the reduced matrix B
cu t (not L 1 AU 1 ). The actual change to the hierarchical basis is defined implicitly
through this recursion. The final hierarchical bases stiffness matrix A 0 is far less sparse,
but generally better conditioned than A, so that standard iterative methods can be
effectively applied to linear systems involving A 0 . In the case of classical HBMG, the
iterative method is just a standard block symmetric Gauss-Seidel iteration, with the
blocks defined in terms of refinement levels in the mesh. See [6, 7] for details.
The multigraph method, ILU-MG, replaces this Gauss-Seidel iteration with an
incomplete factorization. In particular we set
where L is unit lower triangular, D is diagonal, U is unit upper triangular, and E is
the so-called error matrix. The sparsity pattern of L+D+U is defined in terms of the
elimination algorithm described in section 2. This is precisely the sparsity
pattern recursively generated by the hierarchical transformations defined above. By
using a standard ILU factorization, we avoid the recursion of HBMG. At the same
time, we hope that by allowing this additional fill-in, ILU-MG will inherit the desirable
properties of HBMG as a preconditioner.
1 Here we are speaking in generic terms and in particular are not taking into account the possible
occurrence of so-called accidental zeros.
4. Implementation. Our multigraph algorithm is divided into four distinct
phases, analogous to classical sparse Gaussian elimination algorithms.
1. Ordering: A permutation matrix P is computed to reorder the matrix, PAP t .
In addition, vertex parents for each eliminated vertex are defined and determine
the fill-in pattern for the second phase.
2. Symbolic factorization: The (incomplete) fill-in is computed using the graph
of PAP t and the vertex parents. The output is a static data structure for
the incompletely factored sparse matrix.
3. Numerical factorization: The numerical values of the L, D, and U factors are
computed using a MILU factorization.
4. Solution: The solution of computed using a conjugate or biconjugate
gradient algorithm, preconditioned by the incomplete factorization.
This section is divided into two parts. First, we describe the ordering strategy used
to compute P and the vertex parents for each vertex. Next, our MILU factorization
is discussed. We do not believe that our present algorithms for these two phases are
optimal in any sense. However, they are the best ones we have found so far and their
performance seems to justify further work in this area.
4.1. Ordering. In the case of classical sparse Gaussian elimination, ordering
consists of finding a permutation matrix P such that the reordered matrix PAP t
has some desired property in terms of the ensuing factorization. Normally, the permutation
matrix P is constructed based solely on the graph of the matrix (e.g., a
minimum degree ordering [16, 25, 18]) and not on the values of the matrix elements.
In the multigraph algorithm, both the graph and the numerical values of the matrix
A are used to construct both the ordering and the vertex parents.
To simplify notation, we will describe only how the first vertex is ordered and
its parents are selected. The remaining vertices are ordered by the same algorithm
applied inductively. Let
For
fl is given by
ae
is the number of fill-in edges which must be added if v j 2 adj(v i ) is chosen
as the only vertex parent of v i . For the case of one vertex parent, the quality function
defined by
The tentative vertex parent is the vertex
are broken arbitrarily. The quality function q 1 (v i ) represents a compromise between
choosing a parent v j which is as strongly connected to v i as possible (measured in
terms of the size of the off-diagonal matrix elements), and choosing v j to cause as
little fill-in as possible in terms of the factorization.
The size of the parameter fl can be used to indirectly control the number of
fill-in edges resulting from the ordering. Smaller values for fl result in less fill-in.
Experimentally, we determined to be a good choice.
The quality function q 2 (v i ) for the case of two vertex parents is developed in a
similar fashion. Suppose that v j). Then the qualities s ijk
are
given by
s ijk
ae
Let g ijk denote the number of fill-in edges required if fv are chosen as parents;
then
s ijk
and set
Our two parent algorithm actually offers the possibility to each vertex of having zero,
one, or two parents. The choice of parents is based on maximizing the function q 2 (v i ).
As in the single parent case, the quality function q 2 (v i ) seeks a compromise between
choosing strongly connected parents, and choosing parents which allow low fill-in;
experimentally, we determined for the two parent algorithm.
We now describe our algorithm for ordering vertices and computing vertex par-
ents. We begin by computing the quality q p (v i ) for each vertex in the mesh, using (5)
for the case parent algorithm) or (7) for the case
parent algorithm). Along with the quality function, tentative parents are assigned to
each vertex. In cases where no parents can be assigned, q p (v The vertices are
then placed in a heap with the vertex of highest quality at the root. This vertex,
say v i , is ordered first, and its tentative parents become its actual parents. We then
update the graph and compute the reduced matrix. The quality function and heap
position for vertices are updated, as vertices in this set are the only ones
potentially affected by the elimination of v i . This process continues inductively until
all the vertices are ordered (the usual case), or all remaining vertices have q p (v k
(in which case an arbitrary order is assigned to the remaining vertices). We summarize
this algorithm below:
1. Initialize by computing q p (v i ) and tentative parents for
all vertices in a heap according to q p (v i ).
2. While the heap is not empty and the root vertex has q p (v
below.
i. Order the vertex v i at the root of the heap; update the heap. The
tentative parents of v i become the actual parents.
ii. Eliminate v i from the current graph; add fill-in edges as required.
Update the partially factored matrix using the MILU decomposition.
iii. For update the position of v j in the
heap.
We note that step ii above essentially requires an incomplete factorization of the
matrix A to occur concurrently with the ordering, since the quality function q p is
updated based on the current state of the factorization. This makes the ordering
algorithm rather expensive, often as expensive as or even more expensive than the
actual solution. This is partly because this factorization must use a dynamic data
structure, rather than the static sparse matrix data structures which we employ else-
where. On the other hand, as in general sparse matrix calculations, in many cases
ordering can be done once and then used for several factorizations and solutions.
There are several interesting variations of our ordering algorithm which merit
some discussion. Both are related to step iii above and provide a means of partitioning
the matrix in order to formulate block iterative methods. First, in step iii, we can
(artificially) set q p (v parent of v i . This forces v j (and all
other vertices chosen as actual parents) to reside near the bottom of the heap. When
the heap contains only vertices with q the remaining vertices are called coarse
graph vertices and those eliminated are fine graph vertices. This effectively provides
a two level blocking in a fashion quite analogous to the classical two level HBMG
algorithm. If we (correctly) reinitialize q p (v i ) and compute tentative parents for all
vertices remaining in the heap, and then restart the elimination process, we are led
to a natural multilevel blocking, which could form the basis of a block incomplete
factorization algorithm.
Second, in step iii we can artificially set q p (v
the fine graph (eliminated) vertices will form an independent set, in that the diagonal
block of both the original and factored matrices corresponding to this set of vertices
will be diagonal. Reinitializing and restarting the elimination process as above would
result in some multicolor-like ordering, which might have some interesting applications
for vector or parallel processing.
For either of these alternatives, using an enhanced quality function that includes
additional information about the blocking strategy (e.g., bias q p to favor producing
the largest number of vertices in the fine graph set within the other constraints) seems
appropriate.
4.2. Numerical factorization. Our implementation of the MILU factorization
is defined as follows. Let
and set actually compute an incomplete factorization of the matrix
Making such an a priori shift is one simple way to ensure the existence and
stability of our factorization. To describe the first step of our factorization procedure,
let
c=d I
I
The sparsity pattern of matrix \Gammacr t =d generally will not coincide with the sparsity
pattern we choose to allow. Thus we set \Gammacr t =d has the required
sparsity pattern and E 1 is the error matrix for the first step. If we (inductively)
continue the factorization as B
E, we have
c=d I
I
U
Then we have
The matrix \Gammacr t =d is decomposed as S 1 +E 1 by a procedure similar to classical
MILU [22]. Suppose that p ? q and the element c p r q =d is not allowed in the fill-in
pattern. Consider the matrix F pq which is zero except for the four elements
F pq
We then set
where the sum is taken over all (p; q) pairs falling outside the allowed fill-in pat-
tern. This is a typical MILU approach, although the averaging of diagonal entries
appearing in F pq for the nonsymmetric case is a bit unusual.
5. Numerical results. In this section we present a few numerical results for a
first version of our multigraph algorithm. We compare both the one parent and two
parent versions of multigraph with the well-known algebraic multigrid code AMG of
Ruge, St-uben, and Hempel. A suite of six text problems were constructed using the
PLTMG package, version 7.8 [5].
For each problem, three nonuniform adaptive grids with
were generated. For each test case, both the sparse matrix and the right hand side
were saved in a file to serve as input for the iterative solvers. 2 The specific definition
of each test problem is described 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,
shown in Figure 5, is generally very smooth but has some boundary singularities.
Problem Hole. This problem features discontinuous, anisotropic coefficients. The
domain consists of three subregions. On the inner region, the problem is
. 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, shown in Figure 5,
is also relatively smooth, singularities exist at the internal interfaces.
These files are available upon request.
Problem Texas. This is an indefinite 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 \Theta 16, so this
problem is fairly indefinite, as illustrated in Figure 6.
Problem UCSD. This is a simple constant coefficient convection-diffusion equation
\Gammar
posed on a domain in the shape of the UCSD logo. Homogeneous
Dirichlet boundary conditions are imposed. As seen in Figure 6, boundary layers are
formed at the bottom of the region and the top of the 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-diffusion equation of the form
\Gammar
The domain has seven subregions; in the upper left and large lower region. In the
narrow curved band, and is directed radially. Dirichlet boundary conditions
are imposed along the bottom boundary and along a short
segment on the upper left boundary, respectively. Homogeneous Neumann boundary
conditions are specified elsewhere. The solutions, see Figure 7, vary exponentially
across the domain which is typical of semiconductor problems.
In the first 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.
All problems were run on a SGI R10000 Octane with 256 Mb of memory and
compiled with the Fortran f77 -O -64 options. Each run consisted of two phases. The
setup phase consisted of performing several initialization steps. For multigraph al-
gorithms, this included ordering, symbolic factorization, and numerical factorization.
Of these three steps, the sparse ordering is by far the most dominant. The initialization
step for AMG consisted of determining the multigrid levels and constructing
the interpolation operators as well as the coarse grid matrices. In the second phase,
each problem was solved to a relative accuracy in the residual of 10 \Gamma6 starting from
an initial guess of zero.
The results of this comparison are shown in Table 1 and Table 2. Since there
was little variation in the timings for the setup phase, these times (in seconds) are
averaged over all problems with the same grid size.

Table
Average setup time vs. problem size.
N AMG One parent Two parents
Fig. 5. Lake Superior and Hole problems.
Fig. 6. Indefinite and boundary layer problems.
Fig. 7. Semiconductor convection-diffusion problems.
In

Table

2, we present the number of AMG cycles and multigraph solves. The
multigraph algorithm is used as a preconditioner for the composite step conjugate
gradient (CSCG) procedure for symmetric problems of the composite step biconjugate
gradient (CSBCG) procedure for nonsymmetric problems. We count solves rather
than iterations, since composite steps cost about twice as much as single steps. Also,
solves for nonsymmetric problems are twice as expensive as for symmetric problems,
since both A and A t must be preconditioned. The digits columns refer to
where r k is the residual at the kth iteration (cycle). To compare these methods, we
have chosen the CPU time, measured in seconds, to solve each problem.

Table
Performance comparison.
AMG One parent Two parents
Cycles Digits Time Solves Digits Time Solves Digits Time
Superior
Hole
26 6.3 7.29
Texas
UCSD
Jcn
The timings for the one parent and the two parent versions exhibit a greater
than linear growth as a function of problem size. For conventional test problems on
uniform 5-point grids, we observe an O(N log N) dependence, which is significantly
better than most other ILU methods. In contrast, AMG is observed to be O(N) for
this class of problems. However, for all the problems we considered in this study, our
algorithm is competitive with AMG.

Acknowledgments

. The authors wish to thank Steve McCormick for helpful
comments and John Ruge for his invaluable assistance in running the AMG code.



--R

The algebraic multilevel iteration methods - theory and applications

Algebraic multilevel preconditioning methods I
A class of hybrid algebraic multilevel preconditioning methods
PLTMG: A Software Package for Solving Elliptic Partial Differential Equations

The hierarchical basis multigrid method
General sparse elimination requires no permanent integer storage

The hierarchical basis multigrid method and incomplete LU decomposi- tion

Towards algebraic multigrid for elliptic problems of second order


Black box multigrid
Algorithms and data structures for sparse symmetric Gaussian elimination
Algebraic analysis of the hierarchical basis preconditioner
Computer Solution of Large Sparse Positive Definite Systems
Interpolation and related coarsening techniques for the algebraic multigrid method
An algebraic hierarchical basis preconditioner

Incomplete Decompositions - Theory
Analysis of multilevel graph partitioning

A graph theoretic study of the numeric solution of sparse positive definite systems

On the multi-level splitting of finite element spaces
--TR

--CTR
P. K. Jimack, Domain decomposition preconditioning for parallel PDE software, Engineering computational technology, Civil-Comp press, Edinburgh, UK, 2002
Michele Benzi, Preconditioning techniques for large linear systems: a survey, Journal of Computational Physics, v.182 n.2, p.418-477, November 2002
Randolph E. Bank, Compatible coarsening in the multigraph algorithm, Advances in Engineering Software, v.38 n.5, p.287-294, May, 2007
algebraic multigrid;incomplete LU factorization;hierarchical basis