An Approximate Minimum Degree Ordering Algorithm.
An approximate minimum degree (AMD) ordering algorithm  for preordering a  symmetric sparse matrix  prior to numerical factorization is presented.  We use techniques  based on the quotient graph for matrix factorization that allow us to  obtain computationally cheap bounds for the minimum degree. We show that  these bounds are often equal to the actual degree.  The resulting  algorithm is typically much faster than previous minimum degree ordering  algorithms and produces results that are comparable in quality with the  best orderings from other minimum degree algorithms.
Introduction
. When solving large sparse symmetric linear systems of the
is common to precede the numerical factorization by a symmetric
reordering. This reordering is chosen so that pivoting down the diagonal in order on
the resulting permuted matrix PAP produces much less fill-in and work than
computing the factors of A by pivoting down the diagonal in the original order. This
reordering is computed using only information on the matrix structure without taking
account of numerical values and so may not be stable for general matrices. However,
if the matrix A is positive-definite [21], a Cholesky factorization can safely be used.
This technique of preceding the numerical factorization with a symbolic analysis can
also be extended to unsymmetric systems although the numerical factorization phase
must allow for subsequent numerical pivoting [1, 2, 16]. The goal of the preordering
is to find a permutation matrix P so that the subsequent factorization has the least
fill-in. Unfortunately, this problem is NP-complete [31], so heuristics are used.
The minimum degree ordering algorithm is one of the most widely used heuristics,
since it produces factors with relatively low fill-in on a wide range of matrices. Because
of this, the algorithm has received much attention over the past three decades. The
algorithm is a symmetric analogue of Markowitz' method [26] and was first proposed
by Tinney and Walker [30] as algorithm S2. Rose [27, 28] developed a graph theoretical
model of Tinney and Walker's algorithm and renamed it the minimum degree
algorithm, since it performs its pivot selection by selecting from a graph a node of
minimum degree. Later implementations have dramatically improved the time and
memory requirements of Tinney and Walker's method, while maintaining the basic
idea of selecting a node or set of nodes of minimum degree. These improvements have
reduced the memory complexity so that the algorithm can operate within the storage
of the original matrix, and have reduced the amount of work needed to keep track of
the degrees of nodes in the graph (which is the most computationally intensive part
Computer and Information Sciences Department University of Florida, Gainesville, Florida,
USA. phone: (904) 392-1481, email: davis@cis.ufl.edu. Support for this project was provided by
the National Science Foundation (ASC-9111263 and DMS-9223088). Portions of this work were
supported by a post-doctoral grant from CERFACS.
y ENSEEIHT-IRIT, Toulouse, France. email: amestoy@enseeiht.fr.
z Rutherford Appleton Laboratory, Chilton, Didcot, Oxon. 0X11 0QX England, and European
Center for Research and Advanced Training in Scientific Computation (CERFACS), Toulouse,
France. email: isd@letterbox.rl.ac.uk.
P. AMESTOY, T. A. DAVIS, , AND I. S. DUFF
of the algorithm). This work includes that of Duff and Reid [10, 13, 14, 15]; George
and McIntyre [23]; Eisenstat, Gursky, Schultz, and Sherman [17, 18]; George and Liu
[19, 20, 21, 22]; and Liu [25]. More recently, several researchers have relaxed this
heuristic by computing upper bounds on the degrees, rather than the exact degrees,
and selecting a node of minimum upper bound on the degree. This work includes
that of Gilbert, Moler, and Schreiber [24], and Davis and Duff [7, 8]. Davis and Duff
use degree bounds in the unsymmetric-pattern multifrontal method (UMFPACK), an
unsymmetric Markowitz-style algorithm. In this paper, we describe an approximate
minimum degree ordering algorithm based on the symmetric analogue of the degree
bounds used in UMFPACK.
Section 2 presents the original minimum degree algorithm of Tinney and Walker
in the context of the graph model of Rose. Section 3 discusses the quotient graph
(or element graph) model and the use of that model to reduce the time taken by
the algorithm. In this context, we present our notation for the quotient graph, and
present a small example matrix and its graphs. We then use the notation to describe
our approximate degree bounds in Section 4. The Approximate Minimum Degree
(AMD) algorithm and its time complexity is presented in Section 5. In Section 6,
we first analyse the performance and accuracy of our approximate degree bounds on
a set of test matrices from a wide range of disciplines. The AMD algorithm is then
compared with other established codes that compute minimum degree orderings.
2. Elimination graphs. The nonzero pattern of a symmetric n-by-n matrix,
A, can be represented by a graph G ng and
only if a ij 6= 0. Since A is symmetric, G 0 is
undirected.
The elimination graph, G describes the nonzero pattern of the sub-matrix
still to be factorized after the first k pivots have been chosen. It is undirected,
since the matrix remains symmetric as it is factorized. At step k, the graph G k depends
on G k\Gamma1 and the selection of the kth pivot. To find G k , the kth pivot node p
is selected from V k\Gamma1 . Edges are added to E k\Gamma1 to make the nodes adjacent to p in
G k\Gamma1 a clique (a fully connected subgraph). This addition of edges (fill-in) means that
we cannot know the storage requirements in advance. The edges added correspond to
fill-in caused by the kth step of factorization. A fill-in is a nonzero entry L ij , where
zero. The pivot node p and its incident edges are then removed from the
graph G k\Gamma1 to yield the graph G k . Let Adj G k(i) denote the set of nodes adjacent to
i in the graph G k . Throughout this paper, we will use the superscript k to denote a
graph, set, or other structure obtained after the first k pivots have been chosen. For
simplicity, we will drop the superscript when the context is clear.
The minimum degree algorithm selects node p as the kth pivot such that the
degree of p, t p j jAdj G k\Gamma1 (p)j, is minimized (where j:::j denotes the size of a set or the
number of nonzeros in a matrix). The minimum degree algorithm is a non-optimal
greedy heuristic for reducing the number of new edges (fill-ins) introduced during the
factorization. We have already noted that the optimal solution is NP-complete [31].
By minimizing the degree, the algorithm minimizes the upper bound on the fill-in
caused by the kth pivot. Selecting p as pivot creates at most (t 2
new edges in
G.
3. Quotient graphs. In contrast to the elimination graph, the quotient graph
models the factorization of A using an amount of storage that never exceeds the
storage for the original graph, G 0 [21]. The quotient graph is also referred to as the
generalized element model [13, 14, 15, 29]. An important component of a quotient
graph is a clique. It is a particularly economic structure since a clique is represented
by a list of its members rather than by a list of all the edges in the clique. Following
the generalized element model, we refer to nodes removed from the elimination graph
as elements (George and Liu refer to them as eliminated nodes). We use the term
variable to refer to uneliminated nodes.
The quotient graph, G
implicitly represents the elimination
graph G k , where G . For clarity, we drop the superscript k in the following. The
nodes in G consist of variables (the set V ), and elements (the set V ). The edges are
divided into two sets: edges between variables E ' V \Theta V , and between variables and
elements . There are no edges between elements since they are removed
by element absorption. The sets V 0 and E 0 are empty.
We use the following set notation (A, E , and L) to describe the quotient graph
model and our approximate degree bounds. Let A i be the set of variables adjacent to
variable i in G, and let E i be the set of elements adjacent to variable i in G (we refer
to E i as element list i). That is, if i is a variable in V , then
and
The set A i refers to a subset of the nonzero entries in row i of the original matrix A
(thus the notation A). That is, A 0
Let L e denote the set of variables adjacent to element e in G. That is, if e is an
element in V , then we define
The edges E and E in the quotient graph are represented explicitly as the sets A i
each variable in G, and the sets L e for each element in G. We will use A,
and L to denote three sets containing all A i , respectively, for all variables
i and all elements e. George and Liu [21] show that the quotient graph takes no more
storage than the original graph (jA
The quotient graph G and the elimination graph G are closely related. If i is a
variable in G, it is also a variable in G, and
where the "n" is the standard set subtraction operator.
When variable p is selected as the kth pivot, element p is formed (variable p is removed
from V and added to V ). The set L using Equation (3.1).
The set L p represents a permuted nonzero pattern of the kth column of L (thus the
notation L). If i 2 L p , where p is the kth pivot, and variable i will become the mth
pivot (for some m ? k), then the entry Lmk will be nonzero. Equation (3.1) implies
that L e n fpg ' L p for all elements e adjacent to variable p. This means that all variables
adjacent to an element e are adjacent to the element p and these elements
4 P. AMESTOY, T. A. DAVIS, , AND I. S. DUFF
Fig. 3.1. Elimination graph, quotient graph, and matrix for first three steps.45810758104581013692136921369273692357108357108357108357108458102
G 1G
(a) Elimination graph
(b) Quotient graph
(c) Factors and active submatrix
are no longer needed. They are absorbed into the new element p and deleted
[15], and reference to them is replaced by reference to the new element p. The new
element p is added to the element lists, E i , for all variables i adjacent to element p.
Absorbed are removed from all element lists. The sets A p and E p ,
and L e for all e in E p , are deleted. Finally, any entry j in A i , where both i and j are
in L p , is redundant and is deleted. The set A i is thus disjoint with any set L e for
In other words, A k
i is the pattern of those entries in row i of A that are not
modified by steps 1 through k of the Cholesky factorization of PAP T . The net result
is that the new graph G takes the same, or less, storage than before the kth pivot was
selected.
3.1. Quotientgraph example. We illustrate the sequence of elimination graphs
and quotient graphs of a 10-by-10 sparse matrix in Figures 3.1 and 3.2. The example
is ordered so that a minimum degree algorithm recommends pivoting down the diagonal
in the natural order (that is, the permutation matrix is the identity). In Figures
3.1 and 3.2, variables and elements are shown as thin-lined and heavy-lined circles,
respectively. In the matrices in these figures, diagonal entries are numbered, unmodified
original nonzero entries (entries in are shown as a solid squares. The solid
squares in row i form the set A i . The variables in current unabsorbed elements (sets
are indicated by solid circles in the columns of L corresponding to the unabsorbed
Fig. 3.2. Elimination graph, quotient graph, and matrix for steps 4 to 7.36925810
G 6
(a) Elimination graph
(c) Factors and active submatrix
(b) Quotient graph
elements. The solid circles in row i form the set E i . Entries that do not correspond
to edges in the quotient graph are shown as an \Theta. Figure 3.1 shows the elimination
graph, quotient graph, and the matrix prior to elimination (in the left column) and
after the first three steps (from left to right). Figure 3.2 continues the example for
the next four steps.
Consider the transformation of the graph G 2 to the graph G 3 . Variable 3 is
selected as pivot. We have simple case of Equation (3.1)). The
new element 3 represents the pairwise adjacency of variables 5, 6, and 7. The explicit
edge (5,7) is now redundant, and is deleted from A 5 and A 7 .
Also consider the transformation of the graph G 4 to the graph G 5 . Variable 5 is
selected as pivot. The set A 5 is empty and 3g. Following Equation (3.1),
which is the pattern of column 5 of L (excluding the diagonal). Since the new element
5 implies that variables 6, 7, and, 9 are pairwise adjacent, elements 2 and 3 do not
add any information to the graph. They are removed, having been "absorbed" into
element 5. Additionally, the edge (7, is redundant, and is removed from A 7 and
6 P. AMESTOY, T. A. DAVIS, , AND I. S. DUFF
A 9 . In G 4 , we have
After these transformations, we have in G 5 ,
and the new element in G 5 ,
3.2. Indistinguishable variables and external degree. Two variables i and
are indistinguishable in G if Adj G (i) [ They will have the
same degree until one is selected as pivot. If i is selected, then j can be selected
next without causing any additional fill-in. Selecting i and j together is called mass
elimination [23]. Variables i and j are replaced in G by a supervariable containing
both i and j, labeled by its principal variable (i, say) [13, 14, 15]. Variables that
are not supervariables are called simple variables. In practice, new supervariables are
constructed at step k only if both i and j are in L p (where p is the pivot selected at
step k). In addition, rather than checking the graph G for indistinguishability, we use
the quotient graph G so that two variables i and j are found to be indistinguishable
if Adj G (i) [ This comparison is faster than determining if
two variables are indistinguishable in G, but may miss some identifications because,
although indistinguishability in G implies indistinguishability in G, the reverse is not
true.
We denote the set of simple variables in the supervariable with principal variable
i as i, and define is a simple variable. When p is selected as pivot at the
kth step, all variables in p are eliminated. The use of supervariables greatly reduces
the number of degree computations performed, which is the most costly part of the
algorithm. Non-principal variables and their incident edges are removed from the
quotient graph data structure when they are detected. The set notation A and L refers
either to a set of supervariables or to the variables represented by the supervariables,
depending on the context. In degree computations and when used in representing
elimination graphs, the sets refer to variables; otherwise they refer to supervariables.
In

Figure

3.2, detected supervariables are circled by dashed lines. Non-principal
variables are left inside the dashed supervariables. These are, however, removed from
the quotient graph. The last quotient graph in Figure 3.2 represents the selection of
pivots 7, 8, and 9, and thus the right column of the figure depicts G 7 , G 9 , and the
matrix after the ninth pivot step.
The external degree d of a principal variable i is
since the set A i is disjoint from any set L e for
fill-ins occur
if all variables in i are selected as pivots. We refer to t i as the true degree of variable i.
Selecting the pivot with minimum external degree tends to produce a better ordering
than selecting the pivot with minimum true degree [25] (also see Section 6.2).
Algorithm 1 (Minimum degree algorithm, based on quotient graph)
to n do
end for
while k - n do
mass elimination:
select variable that minimizes d p
for each i 2 L p do
remove redundant entries:
element absorption:
compute external degree:
end for
supervariable detection, pairs found via hash function:
for each pair i and j 2 L p do
if i and j are indistinguishable then
remove the supervariable j:
end for
convert variable p to element p:
while
3.3. Quotient-graph-based minimum degree algorithm. A minimum degree
algorithm based on the quotient graph is shown in Algorithm 1. It includes
element absorption, mass elimination, supervariables, and external degrees. Super-
variable detection is simplified by computing a hash function on each variable, so that
not all pairs of variables need be compared [3]. Algorithm 1 does not include two
8 P. AMESTOY, T. A. DAVIS, , AND I. S. DUFF
important features of Liu's Multiple Minimum Degree algorithm (MMD): incomplete
update [17, 18] and multiple elimination [25]. With multiple elimination, an independent
set of pivots with minimum degree is selected before any degrees are updated. If
a variable is adjacent to two or more pivot elements, its degree is computed only once.
A variable j is outmatched if Adj G (i) ' Adj G (j). With incomplete degree update,
the degree update of the outmatched variable j is avoided until variable i is selected
as pivot. These two features further reduce the amount of work needed for the degree
computation in MMD. We will discuss their relationship to the AMD algorithm in
the next section.
The time taken to compute d i using Equation (3.2) by a quotient-graph-based
minimum degree algorithm is
which is \Omega\Gamma jAdj G k(i)j) if all variables are simple. 1 This degree computation is the
most costly part of the minimum degree algorithm. When supervariables are present,
the time taken is in the best case proportional to the degree of the variable in the
"compressed" elimination graph, where all non-principal variables and their incident
edges are removed.
4. Approximate degree. Having now discussed the data structures and the
standard minimum degree implementations, we now consider our approximation for
the minimum degree and indicate its lower complexity.
We assume that p is the kth pivot, and that we compute the bounds only for
supervariables than computing the exact external degree, d i , our
Approximate Minimum Degree algorithm (AMD) computes an upper bound [7, 8],
d
The first two terms (n \Gamma k, the size of the active submatrix, and d
the worst case fill-in) are usually not as tight as the third term in Equation (4.1).
Algorithm 2 computes jL e n L p j for all elements e in the entire quotient graph. The
set L e splits into two disjoint subsets: the external subset L e n L p and the internal
subset scans element e, the term w(e) is initialized to jL e j
and then decremented once for each variable i in the internal subset L e " L p , and, at
the end of Algorithm 2, we have
does not scan element e, the term w(e) is less than zero. Combining these two cases,
we obtain
ae w(e) if w(e) - 0
oe
1 Asymptotic complexity notation is defined in [6]. We write there exist
positive constants c 1 , c 2 , and n 0 such that Similarly,
there exist positive constants c and n 0 such that 0 - cg(n) - f(n) for all n ?
and there exist positive constants c and n 0 such that 0 - f(n) - cg(n) for all
Algorithm 2 (Computation of jL e n
assume
for each supervariable i 2 L p do
for each element e do
end for
end for
Algorithm 2 is followed by a second loop to compute our upper bound degree, d i
for each supervariable i 2 L p , using Equations (4.1) and (4.2). The total time for
Algorithm 2 is
The second loop to compute the upper bound degrees takes time
which is thus equal to the total asymptotic time.
Multiple elimination [25] improves the minimumdegree algorithm by updating the
degree of a variable only once for each set of independent pivots. Incomplete degree
update [17, 18] skips the degree update of outmatched variables. We cannot take full
advantage of the incomplete degree update since it avoids the degree update for some
supervariables adjacent to the pivot element. With our technique (Algorithm 2), we
must scan the element lists for all supervariables i in L p . If the degree update of one
of the supervariables is to be skipped, its element list must still be scanned so that the
external subset terms can be computed for the degree update of other supervariables
in L p . The only advantage of multiple elimination or incomplete degree update would
be to skip the second loop that computes the upper bound degrees for outmatched
variables or supervariables for which the degree has already been computed.
If the total time in Equation (4.3) is amortized across the computation of all
supervariables then the time taken to compute d i is
which is \Theta(jAdj G k(i)j) if all variables are simple. Computing our bound takes time
proportional to the degree of the variable in the quotient graph, G. This is much faster
than the time taken to compute the exact external degree (see Equation (3.3)).
4.1. Accuracy of our approximate degrees. Gilbert, Moler, and Schreiber
[24] also use approximate external degrees that they can compute in the same time
as our degree bound d. In our notation, their bound b
d i is
P. AMESTOY, T. A. DAVIS, , AND I. S. DUFF
Since many pivotal variables are adjacent to two or fewer elements when selected,
Ashcraft and Eisenstat [4] have suggested a combination of b
d and d,
e
ae d if jE
Computing e
d takes the same time as d or b
d, except when jE 2. In this case, it
takes O(jA time to compute e
d, whereas computing d or b
d takes \Theta(jA i
In the Yale Sparse Matrix Package [17] the jL e n L p j term for the case
is computed by scanning L e once. It is then used to compute d i for all i 2 L p for
which pg. This technique can also be used to compute e
d, and thus the time
to compute e
d is O(jA
Theorem 1: Relationship between external degree and the three approximate
degree bounds. The equality d
The inequality d
2. Finally, the inequality
2.
Proof:
The bound b
d i is equal to the exact degree when variable i is adjacent to at most
one element (jE 1). The accuracy of their bound is unaffected by the size of A i ,
since entries are removed from A that fall within the pattern L of an element. Thus,
if there is just one element (the current element p, say), the bound b
d i is tight. If jE i j
is two (the current element, p, and a prior element e, say), we have
The bound b
counts entries in the set (L e "
will be an
overestimate in the possible (even likely) case that a variable j exists that is
adjacent to both e and p. Combined with the definition of e
d, we have d
when
2.
If our bound d i is exact for the same reason that b
d i is exact. If jE i j is two
we have
No entry is in both A i and any element L, since these redundant entries are removed
from A i . Any entry in L p does not appear in the external subset (L e n L p ). Thus, no
entry is counted twice, and d 2. Finally, consider both d i and b
when 2. We have
and
Since these degree bounds are only used when computing the degree of a supervariable
2.
Combining the three inequalities in Theorem 1, the inequality d i - d i - e
holds for all values of jE i j.
Note that, if a variable i is adjacent to two elements or less then our bound is
equal to the exact external degree. This is very important, since most variables of
minimum degree are adjacent to two elements or less. Additionally, our degree bounds
take advantage of element absorption, since the bound depends on jE i j after elements
are absorbed.
4.2. Degree computation example. We illustrate the computation of our approximate
external degree bound in Figures 3.1 and 3.2. Variable 6 is adjacent to three
elements in G 3 and G 4 . All other variables are adjacent to two or less elements. In
G 3 , the bound d 6 is tight, since the two sets jL 1 n L 3 j and jL 2 n L 3 j are disjoint.
In graph G 4 , the current pivot element is We compute
5:
The exact external degree of variable 6 is d 4, as can be seen in the elimination
graph G 4 on the left of Figure 3.2(a). Our bound is one more than the exact external
degree, since the variable 5 appears in both L 2 n L 4 and L 3 n L 4 , but is one less than
the bound b
d i which is equal to 6 in this case. Our bound on the degree of variable
6 is again tight after the next pivot step, since elements 2 and 3 are absorbed into
element 5.
5. The approximate minimum degree algorithm. The Approximate Minimum
Degree algorithm is identical to Algorithm 1, except that the external degree, d i ,
is replaced with d i , throughout. The bound on the external degree, d i , is computed
using Algorithm 2 and Equations (4.1) and (4.2). In addition to absorbing elements
in element with an empty external subset (jL e n L absorbed
into element p, even if e is not adjacent to p. This "aggressive" element absorption
improves the degree bounds by reducing jE j.
As in many other minimum degree algorithms, we use linked lists to assist the
search for a variable of minimum degree. List d holds all supervariables i with degree
bound d d. Maintaining this data structure takes time proportional to the total
number of degree computations, or O(jLj).
Computing the pattern of each pivot element, L p , takes a total of O(jLj) time
overall, since each element is used in the computation of at most one other element,
and the total sizes of all elements constructed is O(jLj).
The AMD algorithm is based on the quotient graph data structure used in the
MA27 minimum degree algorithm [13, 14, 15]. Initially, the sets A are stored, followed
by a small amount of elbow room. When the set L p is formed, it is placed in the elbow
room (or in place of A p if jE collection occurs if the elbow room is
exhausted. During garbage collection, the space taken by A i and E i is reduced to
exactly jA for each supervariable i (which is less than or equal to jA 0
the extra space is reclaimed. The space for A e and E e for all elements fully
reclaimed, as is the space for L e of any absorbed elements e. Each garbage collection
P. AMESTOY, T. A. DAVIS, , AND I. S. DUFF
takes time that is proportional to the size of the workspace (normally \Theta(jAj)). In
practice, elbow room of size n is sufficient.
During the computation of our degree bounds, we compute the following hash
function for supervariable detection [3],
which increases the degree computation time by a small constant factor. We place
each supervariable i in a hash bucket according to Hash(i), taking time O(jLj) overall.
If two or more supervariables are placed in the same hash bucket, then each pair of
supervariables i and j in the hash bucket are tested for indistinguishability. If the
hash function results in no collisions then the total time taken by the comparison is
O(jAj).
Ashcraft [3] uses this hash function as a preprocessing step on the entire matrix
(without the mod(n \Gamma 1) term, and with an O(jV j log jV instead of jV j hash
buckets). In contrast, we use this function during the ordering, and only hash those
variables adjacent to the current pivot element.
For example, variables 7, 8, and 9 are indistinguishable in G 5 , in Figure 3.2(a).
The AMD algorithm would not consider variable 8 at step 5, since it is not adjacent
to the pivot element 5 (refer to quotient graph G 5 in Figure 3.2(b)). AMD would
not construct 9g at step 5, since 7 and 9 are distinguishable in G 5 . It would
construct 9g at step 6, however.
The total number of times the approximate degree d i of variable i is computed
during elimination is no more than the number of nonzero entries in row k of L, where
variable i is the kth pivot. The time taken to compute d i is O(jA 0
i j), or equivalently
the number of nonzero entries in row k of the permuted matrix.
The total time taken by the entire AMD algorithm is thus bounded by the degree
computation,
O
This bound assumes no (or few) supervariable hash collisions and a constant number of
garbage collections. In practice these assumptions seem to hold, but the asymptotic
time would be higher if they did not. In many problem domains, the number of
nonzeros per row of A is a constant, independent of n. For matrices in these domains,
our AMD algorithm takes time O(jLj) (with the same assumptions).
6. Performance results. In this section, we present the results of our experiments
with AMD on a wide range of test matrices. We first compare the degree
computations discussed above (t, d, d, e
d, and b
d), as well as an upper bound on the
true degree, t j d We then compare the AMD algorithm with other established
minimum degree codes (MMD and MA27).
6.1. Test Matrices. We tested all degree bounds and codes on all matrices
in the Harwell/Boeing collection of type PUA, RUA, PSA, and RSA [11, 12] (at
orion.cerfacs.fr or numerical.cc.rl.ac.uk), all non-singular matrices in Saad's
collection (at ftp.cs.umn.edu), all matrices in the University of Florida
collection (available from ftp.cis.ufl.edu in the directory pub/umfpack/matrices),

Table
Selected matrices in test set
Matrix n nz Percentage of Description
RAEFSKY3 21,200 733,784 0.00 13.4 fluid/structure interaction, turbulence
VENKAT01 62,424 827,684 0.71 15.7 unstructured 2D Euler solver
BCSSTK32 44,609 985,046 0.20 27.3 structural eng., automobile chassis
developing pipe flow (turbulent)
CT20STIF 52,329 1,323,067 0.77 33.2 structural eng., CT20 engine block
NASASRB 54,870 1,311,227 0.06 35.0 shuttle rocket booster
OLAF 16,146 499,505 0.41 35.2 NASA test problem
RAEFSKY1 3,242 145,517 0.00 38.9 incompressible flow, pressure-driven pipe
CRYSTK03 24,696 863,241 0.00 40.9 structural eng., crystal vibration
RAEFSKY4 19,779 654,416 0.00 41.4 buckling problem for container model
structural eng., crystal vibration
steering mech.
structural eng., automobile component
EX11 16,614 540,167 0.04 43.3 CFD, 3D cylinder & flat plate heat exch.
FINAN512 74,752 261,120 1.32 46.6 economics, portfolio optimization
RIM 22,560 862,411 2.34 63.2 chemical eng., fluid mechanics problem
BBMAT 38,744 1,274,141 5.81 64.4 CFD, 2D airfoil with turbulence
swell problem on square die
chemical eng., light hydrocarbon recovery
chemical eng., light hydrocarbon recovery
ORANI678 2,529 85,426 6.68 86.9 Australian economic model
migration
APPU 14,000 1,789,392 15.64 94.4 NASA test problem (random matrix)
and several other matrices from NASA and Boeing. Of those 378 matrices, we present
results below on those matrices requiring 500 million or more floating-point operations
for the Cholesky factorization, as well as the ORANI678 matrix in the Harwell/Boeing
collection and the EX19 in Saad's collection (a total of 26 matrices). The latter two
are best-case and worst-case examples from the set of smaller matrices.
For the unsymmetric matrices in the test set, we first used the maximum transversal
algorithm MC21 from the Harwell Subroutine Library [9] to reorder the matrix so
that the permuted matrix has a zero-free diagonal. We then formed the symmetric
pattern of the permuted matrix plus its transpose. This is how a minimum degree
ordering algorithm is used in MUPS [1, 2]. For these matrices, Table 6.1 lists the
statistics for the symmetrized pattern.

Table

6.1 lists the matrix name, the order, the number of nonzeros in lower
triangular part, two statistics obtained with an exact minimum degree ordering (using
d), and a description. In column 4, we report the percentage of pivots p such that
2. Column 4 shows that there is only a small percentage of pivots selected
using an exact minimum degree ordering that have more than two elements in their
adjacency list. Therefore, we can expect a good quality ordering with an algorithm
based on our approximate degree bound. In column 5, we indicate how often a
degree d i is computed when jE (as a percentage of the total number of degree
updates). Table 6.1 is sorted according to this degree update percentage. Column 5
thus reports the percentage of "costly" degree updates performed by a minimum
degree algorithm based on the exact degree. For matrices with relatively large values
14 P. AMESTOY, T. A. DAVIS, , AND I. S. DUFF
in column 5, significant time reductions can be expected with an approximate degree
based algorithm.
Since any minimum degree algorithm is sensitive to tie-breaking issues, we randomly
permuted all matrices and their adjacency lists 21 times (except for the random
APPU matrix, which we ran only once). All methods were given the same set of 21
randomized matrices. We also ran each method on the original matrix. On some ma-
trices, the original matrix gives better ordering time and fill-in results for all methods
than the best result obtained with the randomized matrices. The overall comparisons
are not however dependent on whether original or randomized matrices are used. We
thus report only the median ordering time and fill-in obtained for the randomized
matrices.
The APPU matrix is a random matrix used in a NASA benchmark, and is thus
not representative of sparse matrices from real problems. We include it in our test
set as a pathological case that demonstrates how well AMD handles a very irregular
problem. Its factors are about 90% dense. It was not practical to run the APPU
matrix 21 times because the exact degree update algorithms took too much time.
6.2. Comparing the exact and approximate degrees. To make a valid comparison
between degree update methods, we modified our code for the AMD algorithm
so that we could compute the exact external degree (d), our bound (d), Ashcraft and
Eisenstat's bound ( e
d), Gilbert, Moler, and Schreiber's bound ( b
d), the exact true degree
(t), and our upper bound on the true degree (t). The six codes based on d, d, e
d,
d, t, and t (columns 3 to 8 of Table 2) differ only in how they compute the degree.
Since aggressive absorption is more difficult when using some bounds than others,
we switched off aggressive absorption for these six codes. The actual AMD code (in
column 2 of Table 2) uses d with aggressive absorption.

Table

6.2 lists the median number of nonzeros below the diagonal in L (in thou-
sands) for each method. Results 20% higher than the lowest median jLj in the table (or
higher) are underlined. Our upper bound on the true degree (t) and the exact true degree
(t) give nearly identical results. As expected, using minimum degree algorithms
based on external degree noticeably improves the quality of the ordering (compare
columns 3 and 7, or columns 4 and 8). From the inequality d - d - e
d, we
would expect a similar ranking in the quality of ordering produced by these methods.

Table

6.2 confirms this. The bound d and the exact external degree d produce nearly
identical results. Comparing the AMD results and the d column, aggressive absorption
tends to result in slightly lower fill-in, since it reduces jE j and thus improves the
accuracy of our bound. The e
d bound is often accurate enough to produce good re-
sults, but can fail catastrophically for matrices with a high percentage of approximate
pivots (see column 4 in Table 6.1). The less accurate b
d bound produces notably worse
results for many matrices.
Comparing all 378 matrices, the median jLj when using d is never more than 9%
higher than the median fill-in obtained when using the exact external degree, d (with
the exception of the FINAN512 matrix). The fill-in results for d and d are identical
for nearly half of the 378 matrices. The approximate degree bound d thus gives a very
reliable estimation of the degree in the context of a minimum degree algorithm.
The FINAN512 matrix is highly sensitive to tie-breaking variations. Its graph
consists of two types of nodes: "constraint" nodes and "linking" nodes [5]. The
constraint nodes form independent sparse subgraphs, connected together via a tree of
linking nodes. This matrix is a pathological worst-case matrix for any minimumdegree

Table
Median fill-in results of the degree update methods
Matrix Number of nonzeros below diagonal in L, in thousands
AMD d d e
OLAF 2860 2858 2860 2860 3271 3089 3090
RAEFSKY4 7685 7685 7685 7685 9294 8196 8196
BBMAT 19673 19880 19673 21422 37820 21197 21445
APPU 87648 87613 87648 87566 87562 87605 87631
method. All constraint nodes should be ordered first, but linking nodes have low
degree and tend to be selected first, which causes high fill-in. Using a tree dissection
algorithm, Berger, Mulvey, Rothberg, and Vanderbei [5] obtain an ordering with only
1.83 million nonzeros in L.

Table

6.3 lists the median ordering time (in seconds on a SUN SPARCstation
for each method. Ordering time twice that of the minimum median ordering
time listed in the table (or higher) is underlined. Computing the b
d bound is often
the fastest, since it requires a single pass over the element lists instead of the two
passes required for the d bound. It is, however, sometimes slower than d because it
can generate more fill-in, which increases the ordering time (see Equation 5.1). The
ordering time of the two exact degree updates (d and t) increases dramatically as the
percentage of "costly" degree updates increases (those for which jE
Garbage collection has little effect on the ordering time obtained. In the above
runs, we gave each method elbow room of size n. Usually a single garbage collection
occurred. At most two garbage collections occurred for AMD, and at most three for
the other methods (aggressive absorption reduces the memory requirements).
6.3. Comparing algorithms. In this section, we compare AMD with two other
established minimum degree codes: Liu's Multiple Minimum Degree (MMD) code
[25] and Duff and Reid's MA27 code [15]. MMD stores the element patterns L in a
fragmented manner and requires no elbow room [20, 21]. It uses the exact external
degree, d. MMD creates supervariables only when two variables i and j have no
P. AMESTOY, T. A. DAVIS, , AND I. S. DUFF

Table
Median ordering time of the degree update methods
Matrix Ordering time, in seconds
AMD d d e
CT20STIF 6.62 8.66 6.54 7.07 6.31 8.63 6.45
OLAF 1.83 2.56 1.90 2.16 1.83 2.33 1.78
RAEFSKY4 2.32 2.90 2.18 2.45 2.08 3.12 2.07
EX11 2.70 4.06 2.77 3.00 2.60 4.23 2.89
RIM 5.74 10.38 5.69 6.12 5.72 10.01 5.58
APPU 41.75 2970.54 39.83 43.20 40.64 3074.44 38.93
adjacent variables and exactly two adjacent elements
is the current pivot element). A hash function is not required. MMD
takes advantage of multiple elimination and incomplete update.
MA27 uses the true degree, t, and the same data structures as AMD. It detects
supervariables whenever two variables are adjacent to the current pivot element and
have the same structure in the quotient graph (as does AMD). MA27 uses the true
degree as the hash function for supervariable detection, and does aggressive absorp-
tion. Neither AMD nor MA27 take advantage of multiple elimination or incomplete
update.
Structural engineering matrices tend to have many rows of identical nonzero pat-
tern. Ashcraft has found that the total ordering time of MMD can be significantly
improved by detecting these initial supervariables before starting the elimination [3].
We implemented Ashcraft's pre-compression algorithm, and modified MMD to allow
for initial supervariables. We call the resulting code CMMD ("compressed" MMD).
Pre-compression has little effect on AMD, since it finds these supervariables when their
degrees are first updated. The AMD algorithm on compressed matrices together with
the cost of pre-compression was never faster than AMD.

Table

6.4 lists the median number of nonzeros below the diagonal in L (in thou-
sands) for each code. Results 20% higher than the lowest median jLj in the table
(or higher) are underlined. AMD, MMD, and CMMD find orderings of about the
same quality. MA27 is slightly worse because it uses the true degree (t) instead of the
external degree (d).

Table
Median fill-in results of the four codes
Matrix Number of nonzeros below diagonal
in L, in thousands
AMD MMD CMMD MA27
OLAF 2860 2876 2872 3063
BBMAT 19673 19876 19876 21139
APPU 87648 87647 87647 87605
Considering the entire set of 378 matrices, AMD produces a better median fill-in
than MMD, CMMD, and MA27 for 62%, 61%, and 81% of the matrices, respectively.
AMD never generates more than 7%, 7%, and 4% more nonzeros in L than MMD,
CMMD, and MA27, respectively. We have shown empirically that AMD produces an
ordering at least as good as these other three methods for this large test set.
If the apparent slight difference in ordering quality between AMD and MMD is
statistically significant, we conjecture that it has more to do with earlier supervariable
detection (which affects the external degree) than with the differences between the
external degree and our upper bound.

Table

6.5 lists the median ordering time (in seconds on a SUN SPARCstation 10)
for each method. The ordering time for CMMD includes the time taken by the pre-
compression algorithm. Ordering time twice that of the minimum median ordering
time listed in the table (or higher) is underlined. On certain classes of matrices,
typically those from structural analysis applications, CMMD is significantly faster
than MMD. AMD is the fastest method for all but the EX19 matrix. For the other
352 matrices in our full test set, the differences in ordering time between these various
methods is typically less. If we compare the ordering time of AMD with the other
methods on all matrices in our test set requiring at least a tenth of a second of
ordering time, then AMD is slower than MMD, CMMD, and MA27 only for 6, 15,
and 8 matrices respectively. For the full set of matrices, AMD is never more than 30%
slower than these other methods. The best and worst cases for the relative run-time
P. AMESTOY, T. A. DAVIS, , AND I. S. DUFF

Table
Median ordering time of the four codes
Matrix Ordering time, in seconds
AMD MMD CMMD MA27
CT20STIF 6.62 26.00 9.59 9.81
OLAF 1.83 5.67 4.41 2.64
RAEFSKY4 2.32 5.24 2.36 2.91
RIM 5.74 9.09 8.11 10.13
APPU 41.75 5423.23 5339.24 2683.27
of AMD for the smaller matrices are included in Table 6.5 (the EX19 and ORANI678
matrices).
7.

Summary

. We have described a new upper bound for the degree of nodes
in the elimination graph that can be easily computed in the context of a minimum
degree algorithm. We have demonstrated that this upper-bound for the degree is more
accurate than all previously used degree approximations. We have experimentally
shown that we can replace an exact degree update by our approximate degree update
and obtain almost identical fill-in.
An Approximate Minimum Degree (AMD) based on external degree approximation
has been described. We have shown that the AMD algorithm is highly competitive
with other ordering algorithms. It is typically faster than other minimum
degree algorithms, and produces comparable results to MMD (which is also based on
external degree) in terms of fill-in. AMD typically produces better results, in terms
of fill-in and computing time, than the MA27 minimum degree algorithm (based on
true degrees).
8.

Acknowledgments

. We would like to thank John Gilbert for outlining the
portion of the proof to Theorem 1, Joseph Liu for providing a copy of the
MMD algorithm, and Cleve Ashcraft and Stan Eisenstat for their comments on a
draft of this paper.



--R

Factorization of large sparse matrices based on a multifrontal approach
in a multiprocessor environment
Use of level 3 BLAS in the solution of full and
in High Performance Computing:
Symposium on High Performance Computing

Compressed graphs and the minimum degree algorithm
personal communication.

tic programs using tree dissection
Operations Research
Introduction to Algorithms





Journal on Matrix Analysis and Applications in March
On algorithms for obtaining a maximum transversal

Direct Methods for Sparse Matrices
Sparse matrix test problems


A comparison of sparsity orderings for obtaining a pivotal sequence
in Gaussian elimination

Report AERE R10533




Yale sparse matrix
The symmetric codes

Algorithms and data structures
for sparse symmetric Gaussian elimination

A fast implementation of the minimum degree algorithm using
ACM Transactions on Mathematical Software




On the application of the minimum degree algorithm to
finite element systems
Sparse matrices in MATLAB: design and

Modification of the minimum-degree algorithm by multiple elimination



Symmetric Elimination on Sparse Positive Definite Systems and the Potential
PhD thesis


in Graph Theory and Computing
The generalized element method

Direct solutions of sparse network equations by optimally
ordered triangular factorization
Computing the minimum fill-in is NP-complete

Note: all University of Florida technical reports in this list of references are
available in postscript form via anonymous ftp to ftp.
--TR

--CTR
Dany Mezher , Bernard Philippe, Parallel computation of pseudospectra of large sparse matrices, Parallel Computing, v.28 n.2, p.199-221, February 2002
Frank Dellaert , Michael Kaess, Square Root SAM: Simultaneous Localization and Mapping via Square Root                 Information Smoothing, International Journal of Robotics Research, v.25 n.12, p.1181-1203, December  2006
Csaba Mszros, Detecting "dense" columns in interior point methods for linear programs, Computational Optimization and Applications, v.36 n.2-3, p.309-320, April     2007
Abdou Guermouche , Jean-Yves L'excellent, Constructing memory-minimizing schedules for multifrontal methods, ACM Transactions on Mathematical Software (TOMS), v.32 n.1, p.17-32, March 2006
Timothy A. Davis, Algorithm 849: A concise sparse Cholesky factorization package, ACM Transactions on Mathematical Software (TOMS), v.31 n.4, p.587-591, December 2005
Vladimir Rotkin , Sivan Toledo, The design and implementation of a new out-of-core sparse cholesky factorization method, ACM Transactions on Mathematical Software (TOMS), v.30 n.1, p.19-46, March 2004
Timothy A. Davis, A column pre-ordering strategy for the unsymmetric-pattern multifrontal method, ACM Transactions on Mathematical Software (TOMS), v.30 n.2, p.165-195, June 2004
Timothy A. Davis , John R. Gilbert , Stefan I. Larimore , Esmond G. Ng, A column approximate minimum degree ordering algorithm, ACM Transactions on Mathematical Software (TOMS), v.30 n.3, p.353-376, September 2004
Patrick R. Amestoy , Enseeiht-Irit , Timothy A. Davis , Iain S. Duff, Algorithm 837: AMD, an approximate minimum degree ordering algorithm, ACM Transactions on Mathematical Software (TOMS), v.30 n.3, p.381-388, September 2004
Timothy A. Davis , John R. Gilbert , Stefan I. Larimore , Esmond G. Ng, Algorithm 836: COLAMD, a column approximate minimum degree ordering algorithm, ACM Transactions on Mathematical Software (TOMS), v.30 n.3, p.377-380, September 2004
Iain S. Duff, MA57---a code for the solution of sparse symmetric definite and indefinite systems, ACM Transactions on Mathematical Software (TOMS), v.30 n.2, p.118-144, June 2004
Anshul Gupta, Recent advances in direct methods for solving unsymmetric sparse systems of linear equations, ACM Transactions on Mathematical Software (TOMS), v.28 n.3, p.301-324, September 2002
Patrick R. Amestoy , Iain S. Duff , Jean-Yves L'excellent , Xiaoye S. Li, Analysis and comparison of two general sparse solvers for distributed memory computers, ACM Transactions on Mathematical Software (TOMS), v.27 n.4, p.388-421, December 2001
P. Hnon , P. Ramet , J. Roman, PASTIX: a high-performance parallel direct solver for sparse symmetric positive definite systems, Parallel Computing, v.28 n.2, p.301-321, February 2002
Abdou Guermouche , Jean-Yves L'Excellent , Gil Utard, Impact of reordering on the memory of a multifrontal solver, Parallel Computing, v.29 n.9, p.1191-1218, September
Nicholas I. M. Gould , Jennifer A. Scott, A numerical evaluation of HSL packages for the direct solution of large sparse, symmetric linear systems of equations, ACM Transactions on Mathematical Software (TOMS), v.30 n.3, p.300-325, September 2004
Xiaoye S. Li, An overview of SuperLU: Algorithms, implementation, and user interface, ACM Transactions on Mathematical Software (TOMS), v.31 n.3, p.302-325, September 2005
Haifeng Qian , S. S. Sapatnekar, A hybrid linear equation solver and its application in quadratic placement, Proceedings of the 2005 IEEE/ACM International conference on Computer-aided design, p.905-909, November 06-10, 2005, San Jose, CA
Claudson F. Bornstein , Bruce M. Maggs , Gary L. Miller, Tradeoffs between parallelism and fill in nested dissection, Proceedings of the eleventh annual ACM symposium on Parallel algorithms and architectures, p.191-200, June 27-30, 1999, Saint Malo, France
Xiaoye S. Li , James W. Demmel, SuperLU_DIST: A scalable distributed-memory sparse direct solver for unsymmetric linear systems, ACM Transactions on Mathematical Software (TOMS), v.29 n.2, p.110-140, June
Ove Edlund, A software package for sparse orthogonal factorization and updating, ACM Transactions on Mathematical Software (TOMS), v.28 n.4, p.448-482, December 2002
Timothy A. Davis , Iain S. Duff, A combined unifrontal/multifrontal method for unsymmetric sparse matrices, ACM Transactions on Mathematical Software (TOMS), v.25 n.1, p.1-20, March 1999
Nicholas I. M. Gould , Jennifer A. Scott , Yifan Hu, A numerical evaluation of sparse direct solvers for the solution of large sparse symmetric linear systems of equations, ACM Transactions on Mathematical Software (TOMS), v.33 n.2, p.10-es, June 2007
Olaf Schenk , Klaus Grtner, Solving unsymmetric sparse systems of linear equations with PARDISO, Future Generation Computer Systems, v.20 n.3, p.475-487, April 2004
Patrick R. Amestoy , Iain S. Duff , Stphane Pralet , Christof Vmel, Adapting a parallel sparse direct solver to architectures with clusters of SMPs, Parallel Computing, v.29 n.11-12, p.1645-1668, November/December
Michele Benzi, Preconditioning techniques for large linear systems: a survey, Journal of Computational Physics, v.182 n.2, p.418-477, November 2002
ordering algorithms;quotient graph;graph algorithms;approximate minimum degree ordering algorithm;sparse matrices