A column approximate minimum degree ordering algorithm.
Sparse Gaussian elimination with partial pivoting computes the factorization of a sparse matrix A, where the row ordering P is selected during factorization using standard partial pivoting with row interchanges. The goal is to select a column preordering, Q, based solely on the nonzero pattern of A, that limits the worst-case number of nonzeros in the factorization. The fill-in also depends on P, but Q is selected to reduce an upper bound on the fill-in for any subsequent choice of P. The choice of Q can have a dramatic impact on the number of nonzeros in L and U. One scheme for determining a good column ordering for A is to compute a symmetric ordering that reduces fill-in in the Cholesky factorization of ATA. A conventional minimum degree ordering algorithm would require the sparsity structure of ATA to be computed, which can be expensive both in terms of space and time since ATA may be much denser than A. An alternative is to compute Q directly from the sparsity structure of A; this strategy is used by MATLAB's COLMMD preordering algorithm. A new ordering algorithm, COLAMD, is presented. It is based on the same strategy but uses a better ordering heuristic. COLAMD is faster and computes better orderings, with fewer nonzeros in the factors of the matrix.
Introduction
Sparse Gaussian elimination with partial pivoting computes the factorization
for the sparse nonsymmetric matrix A, where P and Q
are permutation matrices, L is a lower triangular matrix, and U is an upper
triangular matrix. Gilbert and Peierls [30] have shown that sparse partial
pivoting can be implemented in time proportional to the number of
point operations required. The method is used by Matlab when solving a
system of linear equations when A is sparse [27]. An improved
implementation is in a sparse matrix package, SuperLU [11]. The solution
process starts by nding a sparsity-preserving permutation Q. Next, the permutation
P is selected during numerical factorization using standard partial
pivoting with row interchanges. The permutation P is selected without regard
to sparsity. Our goal is to compute a sparsity-preserving permutation
Q solely from the pattern of A such that the LU factorization
remains as sparse as possible, regardless of the subsequent choice of P. Our
resulting code has been incorporated in SuperLU and in Matlab Version 6.
Section 2 provides the theoretical background for our algorithm. In Section
3, we describe symbolic LU factorization (with no column interchanges)
as a precursor to the colamd ordering algorithm presented in Section 4. Section
4 also describes the various metrics for selecting columns that we evaluated
in the design of the code. Experimental results for square nonsymmetric
matrices, rectangular matrices, and symmetric matrices are presented in Section
5. Section 6 presents our conclusions, and describes how to obtain the
colamd and symamd codes.
Our notation is as follows. Set subtraction is denoted by the \n" oper-
ator. We use j:::j to denote either the absolute value of a scalar, the number
of nonzero entries in a matrix, or the size of a set. The usage will
be clear in context. The structure of a matrix A, denoted by Struct(A),
is a set that contains the locations of the nonzero entries in A; that is,
0g. Throughout this paper, and in the algo-
rithms, we assume that exact numerical cancellations do not occur. Some
entries in the matrix may happen to become zero during factorization (be-
cause of accidental cancellation); we still refer to these as \nonzero" entries.
Let A be a given nonsymmetric matrix and assume that it is nonsingular. It
follows from Du [14] that the rows (or the columns) of A can be permuted
so that the diagonal entries of the permuted A are all nonzero. We therefore
assume throughout this paper that the given matrix A has been permuted
accordingly so that it has a zero-free diagonal.
Suppose that L and U are the triangular factors obtained when Gaussian
elimination with partial pivoting is applied to A. That is,
some row permutation P. Consider the symmetric positive denite matrix
A T A, which has a Cholesky factorization A T
with L C being lower
triangular. George and Ng [25] showed that if A has a zero-free diagonal,
then the pattern of LC +L T
includes the patterns of L and U, regardless of
the row permutation used in partial pivoting. We summarize the result in
the following theorem.
Theorem 2.1 (George and Ng [25]) Let A be a nonsingular and non-symmetric
matrix that has a zero-free diagonal. Let L C denote the Cholesky
factor of A T A. Let L and U be the LU factors of A obtained by partial
pivoting. Then Struct(U)  Struct(L T
), and the entries within each column
of L can be rearranged to yield a triangular matrix ^
L with Struct(
Gilbert and Ng [28] also showed that the bound on U is tight when A is
a strong Hall matrix.
Theorem 2.2 (Gilbert and Ng [28]) Let A be a nonsingular and non-symmetric
matrix that has a zero-free diagonal. Assume that A is strong
Hall. Let LC denote the Cholesky factor of A T A. Let L and U be the LU factors
of A obtained by partial pivoting. For any choice of (i;
there exists an assignment of numerical values to the nonzero entries of A
such that U ij 6= 0.
It is well known that the sparsity of L C depends drastically on the way in
which the rows and columns of A T A are permuted [23]. Thus, Theorems 2.1
1 A matrix A is strong Hall if every set of k columns of A, 1  k  n 1, contain at
least rows. A strong Hall matrix is irreducible. See Coleman et al. [6] for
details.
and 2.2 suggest a way to reduce the amount of ll in L and U [25, 26].
Given a matrix A, we rst form the pattern of A T A. Then we compute a
symmetric ordering Q of A T A to reduce the amount of ll in LC . Finally,
we apply the permutation Q to the columns of A.
The main disadvantage with the approach above is the cost of forming
A T A. Even if A is sparse, the product A T A may be dense. Consequently,
the time and storage required to form the product may be high. The primary
goal of this paper is to describe ordering algorithms that compute Q from
the pattern of A without explicitly forming A T A.
3 Symbolic LU factorization
Our column ordering technique is based on a symbolic analysis of the conventional
outer-product formulation of Gaussian elimination of the n-by-n matrix
A with row interchanges to maintain numerical stability. Let A A.
For denote the bottom right (n k)-by-(n
submatrix of A (k 1) after Gaussian elimination is performed. We assume
that the rows and columns of A (k) are labeled from k + 1 to n.
At the k-th step of Gaussian elimination, one nds the entry with the
largest magnitude in column k of A (k 1) and swaps rows to place this entry
on the diagonal. Column k of L is then a scaled version of column k of A (k 1) ,
and row k of U is row k of A (k 1) . The outer product of the column k of L
and the row k of U is subtracted from A (k 1) to give A (k) . The result is the
factorization is the permutation matrix determined by
the row interchanges. When A is sparse, the update using the outer product
may turn some of the zero entries into nonzero. These new nonzero entries
are referred to as ll-in. The amount of ll-in that occurs depends on the
order in which the columns of A are eliminated [23]. Our goal, therefore, is to
nd a column ordering Q of the matrix A prior to numerical factorization, so
that the LU factorization of the column-permuted matrix AQ is sparser than
the LU factorization of the original matrix A, regardless of how P is chosen.
Another application of column orderings is sparse QR factorization [33].
In order to control the ll-in that occurs when A is sparse, we need to
know the nonzero patterns of the matrices A (k) , for each k from 0 to n 1.
From this, we can determine a column ordering Q that attempts to keep the
factorization of AQ sparse as the factorization progresses.
We must compute symbolic patterns L, U , and A (k) that account for all
possible row interchanges that may occur during numerical factorization; the
actual structures of L, U, and A (k) during the factorization of PA with any
specic row permutation P will be subsets of these patterns.
Let L k and U k denote the nonzero patterns of column k of L and row k
of U , respectively. For i > k, let R (k)
i denote the set of column indices of
nonzero entries in row i of A (k) ,
R (k)
Similarly, for j > k, let C (k)
denote the set of row indices of nonzero entries
in column j of A (k) ,
Note that these denitions imply
(1)
for all k, i, and j.
We now describe how to compute L k , U k , and A k in order of increasing
k. We begin with A A.
Using the assumption that A has a zero-free diagonal, it is easy to see
that the pattern of the k-th column of L is simply column k of A (k 1) , so we
may take
regardless of how the nonzero numerical values are interchanged. Consider
the possible nonzero pattern of the k-th row of U after partial pivoting. Any
row i can be selected as the pivot row for which a (k 1)
ik
is nonzero. This means
the potential candidate pivot rows correspond to the set L k . The pattern of
the k-th pivot row is bounded by the union of all candidate pivot rows [26],
so we may take
If an arbitrary row r 2 L k is selected as the pivot row, then R (k 1)
for all i 2 L k , so that U k can accommodate any row interchanges due to
numerical partial pivoting.
In the outer product step, multiples of the pivot row are added to each
row i in the pivot column.
, the pattern of each row i in
of the matrix in A (k) after this outer product step must be contained
in the set U k n fkg. Since we do not know which row is the pivot row, we can
bound the pattern of row i of A (k) with
R (k)
Note that all rows i in L k n fkg now have the same pattern, and we can
save time and space by not storing all of them. Any step that modies one
of these rows will modify all of them, and the rows will again have the same
nonzero pattern. Let fkg. At step p, L k n fkg  L p , and so
R (p)
nfkg. Thus, at step k, we can replace all rows i in
with a single super-row R (k)
r
that represents the pattern all
of the rows i in the set L k nfkg. The integer r is an arbitrary place-holder; in
our symbolic factorization algorithm we choose an arbitrary row index from
the set C (k 1)
. We dene [r] as the set of rows (with size j[r]j) represented by
the super-row r. When the symbolic factorization starts, every row belongs
to exactly one set; i.e., frg. To compute the symbolic update, we
need only place the single representative, r, into the set C (k)
, and we can
remove the discarded rows L k n frg. Sherman [43] introduced the concept of
super-columns in the context of sparse Cholesky factorization.
With these simplications, at step k the sets R and C represent a bound
on the pattern of A (k) for LU factorization. They are also a quotient graph
representation [17, 21] of the matrix obtained at the k-th step of the Cholesky
factorization of A T A. Our column preordering method does not require the
patterns of L and U , so these can be discarded. As a result, the algorithm
shown in Figure 1 requires only O(jAj) space. Superscripts on R and C are
dropped for clarity.
The initial storage required is O(jAj). At step k,
because R r is computed as
R r =@ [
R iA n fkg: (2)
The sets R i for all i in C k n frg are then discarded. When the outer product
is formed, each column C j (where j 2 R r ) reduces in size or stays the same
then from (2) this implies 9i 2 C k such that j 2 R i . From
Symbolic LU factorization algorithm
to n do
determine pattern for pivot row and column:
Let r be an arbitrary index in C k
outer product update:
do
end for
end for
end for

Figure

1: Symbolic LU factorization algorithm
which implies that
j. Thus, each step strictly reduces the total storage required
for the pattern of A (k) .
The total time taken is dominated by the symbolic update. Computing
the pivot rows R r for the entire algorithm takes a total of O(jU j) time, since
fkg is the pattern of the o-diagonal part of row k of U , and since
R r is discarded when it is included in a set union for a subsequent pivot
row. Modifying a single column C j in the symbolic update takes no more
than O(jA j j) time, since jC j j  jA j j. Column j is modied at most once
for each row U k that contains column index j. Thus, column j is modied
at most jU j j times during the entire algorithm. Here U j is the upper bound
pattern on column j of the matrix U for any row permutation P due to
partial pivoting. Thus, the total time taken is
where the nonzero pattern of U can accommodate any row permutation P
due to partial pivoting. This is the same time required to compute the sparse
matrix product AU T , and is typically much less than the time required for
numerical factorization, which is equal to the time required to compute the
sparse matrix product L times U [31],
O
Faster methods exist to compute this upper bound on symbolic factorization
with partial pivoting. George and Ng's method [26] for computing the
patterns of L and U takes only O(jLj the column ordering is
xed. The method does not produce the nonzero patterns of A (k) , however,
so it cannot be used as a basis for the colamd algorithm described in the
next section.
4 The colamd ordering algorithm
We now present a column ordering method that is based on the symbolic
LU factorization algorithm presented in the last section and that has the
same time complexity and storage requirements. At step k, we select a pivot
column c to minimize some metric on all of the candidate pivot columns as
an attempt to reduce ll-in. Columns c and k are then exchanged.
The presence of dense (or nearly dense) rows or columns in A can greatly
aect both the ordering quality and run time. A single dense row in A renders
all our bounds useless; the bound on the nonzero pattern of A (1) is a
completely dense matrix. We can hope that this dense row will not be selected
as the rst pivot row during numerical factorization, and withhold the
row from the ordering algorithm. Dense columns do not aect the ordering
quality, but they do increase the ordering time. A single dense column increases
the ordering time by O(n 2 ). Dense columns are withheld from the
symbolic factorization and ordering algorithm, and placed last in the column
ordering Q. Determining how dense a row or column should be for it to be
ignored is problem dependent. We used the same default threshold used by
Matlab, 50%, which is probably too high for most matrices.
Taking advantage of super-columns can greatly reduce the ordering time.
In the symbolic update step, we look at all columns j in the set R r . If any two
or more columns have the same pattern (tested via a hash function [2]), they
are merged into a single super-column. This saves both time and storage in
subsequent steps. Selecting a super-column c allows us to mass-eliminate [24]
all columns [c] represented by the super-column c, and we skip ahead to step
of symbolic factorization. If the pattern of column or super-column
becomes the same as the pivot column pattern (C after the
update, it can be eliminated immediately. Since the columns [c]
represented by a super-column keep the same nonzero pattern until they are
eliminated, we only need to maintain the degree, or other column selection
metric, for the representative column c.
With super-rows, the size of the set C j diers from the sum of the numbers
of rows represented by the super-rows i in the set C j . Similarly, jR i j is the
number of super-columns in row i, which is not the same as the number of
columns represented by the super-columns in R i . Let kR i k denote the sum
of the numbers of rows represented by the super-columns j in R i ,
Similarly, we dene kC j k as
are the same. Similarly, jC
when there are no super-rows.
By necessity, the choice of the pivot column c is a heuristic, since obtaining
an ordering with minimum ll-in is an NP-complete problem [45]. Several
strategies are possible. Each selects a column c that minimizes some metric,
described below. For most of these methods, we need to compute the initial
metric for each column prior to symbolic factorization and ordering, and then
recompute the metric for each column j in the pivot row pattern R r at step
k.
1. Exact external row degree:
Select c to minimize the size of the resulting pivot row,
R iA n fcg

(We exclude the pivot column itself, thus computing an external row
degree [37].) The exact degree is expensive to compute (its time can
dominate the symbolic factorization time), and our experience with
symmetric orderings is that exact degrees do not provide better orderings
than good approximate degrees (such as those from the ordering
methods amdbar and mc47bd) [1]. We thus did not test this method.
2. Matlab approximate external row degree:
The colmmd column ordering algorithm in Matlab is based on the symbolic
LU factorization algorithm presented in Section 3 [27]. Colmmd
selects as pivot the column c that minimizes a loose upper bound on
the external row degree,
Note that c 2 R i , so kR i n j[c]j. Using this bound for
the symbolic update does not increase the asymptotic time complexity
of the ordering algorithm above that of the symbolic LU factorization
algorithm. Computing the Matlab metric for the initial columns of A
is very fast, taking only O(jAj) time.
3. AMD approximate external row degree:
The AMD algorithm for ordering symmetric matrices prior to a Cholesky
factorization is based on a bound on the external row degree that is
tighter than the Matlab bound [1]. It was rst used in the nonsymmetric-
pattern multifrontal method (UMFPACK), to compute the ordering
during numerical factorization [8, 9]. In the context of a column ordering
algorithm to bound the ll-in for the sparse partial pivoting
method, the bound on kR r k is
R s is the most recent pivot row that modied C c in the symbolic
update, and thus s 2 C c . To compute the AMD metric on the initial
columns of A, we select an arbitrary s 2 C c to compute the bound for
column c. Computing the initial AMD metric takes the same amount of
time as it takes to compute the sparse matrix product AA T . Although
this is costly, it does not increase the asymptotic time of our algorithm.
The time to compute this bound during the symbolic update phase is
asymptotically the same as computing the Matlab bound.
4. Size of the Householder update:
The symbolic LU factorization also computes the pattern of R for a QR
factorization of A [6, 32]. At step k, the size of the Householder update
is kC k k-by-kR r k. The term kC k k is known exactly; the kR r k term can
be computed or approximated using any of the above methods. The
method we tested selected c to minimize the product of kC c k and the
AMD approximation of kR r k in (3). We tested and then discarded this
method, since it gave much worse orderings than the other methods.
5. Approximate Markowitz criterion:
The Markowitz heuristic selects as pivot the entry a ij that minimizes
the product of the degrees of row i and column j, which is the amount of
work required in the subsequent outer product step [38]. The criterion
assumes that the rst k 1 pivot rows and columns have been selected.
In our case, we do not know the exact pivot row ordering yet, so we do
not know the true row or column degrees. In our tests, we selected c to
minimize kC c k times max i2Cc kR i k. This is a tighter upper bound on the
actual pivot row degree if column c is selected as the k-th pivot column
and row i is selected as the k-th pivot row. Although R i does not
bound the pattern R r of the k-th pivot row for arbitrary row partial
does bound the size of the actual pivot row once the
row ordering is selected. Since kC c k and kR i k are known exactly, no
approximations need to be used. We tested and then discarded this
method, since it gave much worse orderings than the other methods.
6. Approximate deciency:
The deciency of a pivot is the number of new nonzero entries that
would be created if that pivot is selected; selecting the pivot with least
deciency leads to a minimum deciency ordering algorithm. Exact
deciency is very costly to compute. Approximate minimum deciency
ordering algorithms have been successfully used for symmetric matrices,
in the context of Cholesky factorization [40, 42]. In an nonsymmetric
context, the deciency of column c can be bounded by
Any new nonzeros are limited to the kC c k-by-kR r k Householder up-
date. Each super-row however, is contained in this submatrix
and thus reduces the possible ll-in by the number of nonzeros it repre-
sents. The kC c k and kR i k terms are known exactly; we used the AMD
approximation for kR r k, from (3).
We tested one variant of approximate deciency [35, 36]. The initialization
phase computes the approximate deciency of all columns, and
the approximate deciency is recomputed any time a column is mod-
ied. No aggressive row absorption is used (discussed below), since
it worsens the approximation. The experimental results were mixed.
Compared to our nal colamd variant, approximate deciency was better
for about 60% of the matrices in our large test suite. When it was
worse than colamd it was sometimes much worse, and vice versa.
As a by-product of the AMD row degree computation, we compute the
sizes of the set dierences kR i n R r k when super-row r is the bound on the
pivot row at step k, for all rows i remaining in the matrix. If we nd that this
set dierence is zero, R i is a subset of row R r . The presence of row i does
not aect the exact row degree, although it does aect the Matlab and AMD
approximations. Row i can be deleted, and absorbed into [r]. We refer to
the deletion of R i when
R r as aggressive row absorption. It is similar to
element absorption, introduced by Du and Reid [17]. If AMD row degrees
are computed as the initial metric, then initial aggressive row absorption can
occur in that phase as well. Aggressive row absorption reduces the run time,
since it costs almost nothing to detect this condition when the AMD metric
is used, and results in fewer rows to consider in subsequent steps.
We tested sixteen variants of colamd, based on all possible combinations
of the following design decisions:
1. Initial metric: Matlab or AMD approximation.
2. Metric computed during the symbolic outer product update: Matlab
or AMD approximation.
3. With or without initial aggressive row absorption.
4. With or without aggressive row absorption during the symbolic factorization

Note that some of the combinations are somewhat costly, since aggressive
row absorption is easy when using the AMD row degrees, but dicult when
using the Matlab approximation.
Since the AMD metric was shown to be superior to the Matlab approximation
in the context of minimum degree orderings for sparse Cholesky
factorization [1], we expected the four variants based solely on the AMD
metric to be superior, so much so that we did not intend to test all sixteen
methods. To our surprise, we found that better orderings were obtained with
an initial Matlab metric and an AMD metric during the symbolic update.
We discovered this by accident. A bug in our initial computation of the AMD
metric resulted in the Matlab approximation being computed instead. This
\bug" gave us better orderings, so we kept it. Using an initial Matlab metric
gave, on average, orderings with 8% fewer
oating-point operations in the
subsequent numerical factorization than using an initial AMD metric. The
initial Matlab metric is also faster to compute.
The version of our ordering algorithm that we recommend, which we now
simply call colamd, uses the initial Matlab metric, the AMD metric for the
symbolic update, no initial aggressive row absorption, and aggressive row absorption
during the symbolic update. Since the AMD metric is incompatible
with both incomplete degree update [18, 19] and multiple elimination [37],
it uses neither strategy. In multiple elimination, a set of pivotal columns is
selected such that the symbolic update of any column in the set does not
aect the other columns in the set. The selection of this set reduces the
number of symbolic updates that must be performed.
Matlab's colmmd ordering algorithm uses the Matlab metric throughout,
multiple elimination with a relaxed threshold, super-rows and super-columns,
mass elimination, and aggressive row absorption. Aggressive row absorption
is not free, as it is in colamd. Instead, an explicit test is made every three
stages of multiple elimination. Super-columns are searched for every three
stages, by default. Rows more than 50% dense are ignored. Since the Matlab
metric can lead to lower quality orderings than an exact metric, options are
provided for selecting an exact external row degree metric, modifying the
multiple elimination threshold, changing the frequency of super-column detection
and aggressive row absorption, and changing the dense row threshold.
Selecting spparms ('tight') in Matlab improves the ordering computed by
colmmd, but at high cost in ordering time.
5 Experimental results
We tested our colamd ordering algorithm with three sets of matrices: square
nonsymmetric matrices, rectangular matrices, and symmetric positive de-
nite matrices. Our test set is the entire University of Florida sparse matrix
collection [7], which includes the Harwell/Boeing test set [15, 16], the linear
programming problems in Netlib at http://www.netlib.org [13], as well as
many other matrices. We exclude complex matrices, nonsymmetric matrices
for which only the pattern was provided, and unassembled nite element ma-
trices. Some matrices include explicit zero entries in the description of their
pattern. Since we ignore numerical cancellation, we included these entries
when nding the ordering and determining the resulting factorization. Each
method is compared by ordering time and quality (number of nonzeros in
the factors, and number of
oating-point operations to compute the factor-
ization). Although we tested all matrices in the collection, we present here
a summary of only some of the larger problems (those for which the best
ordering resulted in a factorization requiring 10 7 operations or more to com-
pute). Complete results are presented in Larimore's thesis [36], available as a
technical report at http://www.cise.ufl.edu/. We performed these experiments
on a Sun Ultra Enterprise 4000/5000 with 2 GB of main memory and
eight 248 Mhz UltraSparc-II processors (only one processor was used). The
code was written in ANSI/ISO C and compiled using the Solaris C compiler
(via Mathwork's mex shell script for interfacing Matlab to software written
in C), with strict ANSI/ISO compliance.
5.1 Square nonsymmetric matrices
For square nonsymmetric matrices, we used colamd to nd a column preordering
Q for sparse partial pivoting. These results are compared with
colmmd with default option settings, and with amdbar [1] applied to the pattern
of A T A (ignoring numerical cancellation, and withholding the same
dense rows and columns from A that were ignored by colamd and colmmd).
The amdbar routine is the same as mc47bd in the Harwell Subroutine Library,
except that it does not perform aggressive row absorption. Both amdbar and
mc47bd use the AMD-style approximate degree. For sparse partial pivoting
with the matrices in our test set, amdbar provides slightly better orderings
than mc47bd. After nding a column ordering, we factorized the matrix AQ
with the SuperLU package [11]. This package was chosen since colamd was
written to replace the column ordering in the current version of SuperLU.
SuperLU is based on the BLAS [12]. 2
We excluded matrices that can be permuted to upper block triangular
form [14] with a substantial improvement in factorization time, for two rea-
sons: (1) our bounds are not tight if the matrix A is not strong Hall, and (2)
SuperLU does not take advantage of reducibility to block upper triangular
form. Such matrices should be factorized by ordering and factorizing each
irreducible diagonal submatrix, after permuting the matrix to block triangular
form. Our resulting test set had 106 matrices, all requiring more than
more operations to factorize.
A representative sample is shown in in Table 1, sorted according to colamd
versus colmmd ordering quality. 3 Table 2 reports the ordering time of colamd,
We used the BLAS routines provided in SuperLU because factorization time is a
secondary measure, and because we had diculty using the Sun-optimized BLAS from a
Matlab-callable driver.
3 To obtain the sample, we sorted the 106 matrices according to the relative
point operation count for the colamd ordering and the colmmd ordering, and then selected
every 8th matrix. We also included twotone, the matrix with the largest dimension in
the test set.

Table

1: Unsymmetric matrices
matrix n nnz description
goodwin 7320 324784 nite-element,
uid dynamics, Navier-Stokes
and elliptic mesh (R. Goodwin)
raefsky2 3242 294276 incompressible
ow in pressure-driven pipe
(A. Raefsky)
twotone 120750 1224224 frequency-domain analysis of nonlinear
analog circuit [20]
computational
uid dynamics (CFD),
Charleston Harbor (Steve Bova)
lhr17c 17576 381975 light hydrocarbon recovery problem [47]
exchanger with
ow
redistribution (D. Averous)
rdist3a 2398 61896 chemical process separation [46]
garon2 13535 390607 2D nite-element, Navier-Stokes (A. Garon)
two-phase
uid
ow (D. Graham)
cavity25 4562 138187 driven cavity, nite-element (A. Chapman)
ex20 2203 69981 2D attenuation of a surface disturbance,
nonlinear CFD (Y. Saad)
wang1 2903 19093 electron continuity, 3D diode [39]
ex14 3251 66775 2D isothermal seepage
ow (Y. Saad)
ex40 7740 458012 3D die swell problem on a square die,
computational
uid dynamics (Y. Saad)
colmmd, and amdbar (the time to compute the pattern of A T A is included
in the amdbar time). For comparison, the last column reports the SuperLU
factorization time, using the colamd ordering. Table 3 gives the resulting
number of nonzeros in L + U, and the
oating-point operations required to
factorize the permuted matrix, for each of the ordering methods. The median
results reported in this paper are for just the representative samples, not the
full test sets. The results for the samples and the full test sets are similar.
Colamd is superior to the other methods for these matrices. Colamd was
typically 3.9 times faster than colmmd and 2.1 times faster than amdbar. For
two matrices, colmmd took more time to order the matrix than SuperLU took
to factorize it. The orderings found by colmmd result in a median increase of

Table

2: Ordering and factorization time, in seconds
matrix colamd colmmd amdbar SuperLU
goodwin 0.31 0.42 1.51 44.65
raefsky2 1.11 1.58 2.11 62.36
twotone 5.88 70.82 14.15 306.85
lhr17c 1.62 15.22 3.08 5.80
ex40 2.34 13.58 2.61 27.21
Median time
relative to - 3.9 2.1 37.1

Table

3: Ordering quality, as factorized by SuperLU
Nonzeros in L +U
matrix colamd colmmd amdbar colamd colmmd amdbar
goodwin 5634 3103 2734 1909 665 498
twotone
ex20 483 589 464 48 86 50
ex14 711 1012 885 91 214 157
ex40 4315 8537 6072 1075 5369 2606
Median result
relative to - 1.10 0.99 - 1.36 1.05
10% in nonzeros in the LU factors and 36% in
oating-point operations, as
compared to colamd. The ordering quality of colamd and amdbar are similar,
although there are large variations in both directions in a small number of
matrices. For a few matrices (in our larger test set, not in Table 1) amdbar
requires more space to store A T A than SuperLU requires to factorize the
permuted matrix.
5.2 Rectangular matrices
For m-by-n rectangular matrices with m > n, we found a column ordering
Q for the Cholesky factorization of (AQ) T (AQ), which is one method for
solving a least squares problem. If m < n, we found a row ordering P for
the Cholesky factorization of (PA)D 2 (PA) T , which arises in interior point
methods for solving linear programming problems [34, 44]. Here, D is a
diagonal matrix. In the latter case, colamd and colmmd found a column
ordering of A T and we used that as the row ordering P. We compared these
two methods with amdbar on the corresponding matrix, A T A (if m > n) or
Our test set was limited. It included only 37 matrices requiring more
but three of these were Netlib linear programming
problems (with m < n). A representative selection 4 is shown in Table 4.
We compared the three methods based on their ordering time, number
of nonzeros in the factor L, and
oating-point operation count required for
the Cholesky factorization. These metrics were obtained from symbfact in
Matlab, a fast symbolic factorization [22, 23, 29]. We did not perform the
numerical Cholesky factorization. The time to construct A T A or AA T is
included in the ordering time for amdbar. The results are shown in Tables
5 and 6.
For these matrices, colamd was twice as fast as colmmd and slightly faster
than amdbar. All three produced comparable orderings. Each method has a
few matrices it does well on, and a few that it does poorly on. Colamd and
colmmd both typically require less storage than amdbar, sometimes by several
orders of magnitude, because of the need to compute the pattern of AA T
prior to ordering it with amdbar. The size of the Cholesky factors always
dominates the size of AA T , however. The primary advantage of colamd and
4 Every 5th matrix from the 37 large matrices, in increasing order of colamd versus
colmmd ordering quality, was chosen.

Table

4: Rectangular matrices (all linear programming problems)
gran 2629 2525 20111 60511 British Petroleum operations
lower bound for a
quadratic assignment problem [41]
tting linear inequalities to
data, min. sum of piecewise-linear
penalties (R. Fourer)
pds 20 33874 108175 232647 320196 military airlift operations [5]
d2q06c 2171 5831 33081 56153 (J. Tomlin)

Table

5: Ordering time, in seconds
matrix colamd colmmd amdbar
gran 0.04 0.07 0.06
pds 20 3.86 10.53 3.92
Median time
relative to - 2.00 1.14
colamd

Table

Ordering quality
Nonzeros in L (10 3 ) Flop. count (10 6 )
matrix colamd colmmd amdbar colamd colmmd amdbar
gran 191 158 143 48 33 27
pds 20 6827 8159 6929 8598 12640 8812
Median result
relative to - 1.00 0.99 - 1.02 0.98
colamd
colmmd over amdbar in this context is the ability to analyze the Cholesky
factorization by nding the ordering and size of the resulting factor L in
space proportional to the number of nonzeros in the matrix A, rather than
proportional to the size of the matrix to be factorized
5.3 Symmetric matrices
For symmetric matrices, Matlab's symmmd routine constructs a matrix M
such that the pattern of M T M is the same as A, and then nds a column
ordering of M using colmmd. There is one row in M for each entry a ij below
the diagonal of A, with nonzero entries in column i and j. This method gives
a reasonable ordering for the Cholesky factorization of A. We implemented
an analogous symamd routine that is based on colamd.
Our test set included 50 matrices requiring 10 7 or more operations to
factorize. The largest was a computational
uid dynamics problem from Ed
Rothberg requiring 136 billion operations to factorize (cfd2). Our sample
test matrices 5 are shown in Table 7. Results from symamd, symmmd, amdbar,
and mc47bd are shown in Tables 8 and 9. The time to construct M is included
in the ordering time for symamd and symmmd.
For these 9 matrices, symamd was over six times faster than symmmd, on
5 Every 7th matrix, in order of increasing ratio of symamd to symmmd ordering quality,
was chosen.

Table

7: Symmetric matrices
matrix n nnz description
pwt 36519 326107 pressurized wind tunnel (R. Grimes)
msc00726 726 34518 symmetric test matrix from MSC/NASTRAN
bcsstk38 8032 355460 stiness matrix, airplane engine component (R. Grimes)
engineering problem (NASA Langley)
cfd2 123440 3087898 computational
uid dynamics (E. Rothberg)
3dtube pressure tube (E. Rothberg)
gearbox 153746 9080404 ZF aircraft
ap actuator (E. Rothberg)
finan512 74752 596992 portfolio optimization [4]

Table

8: Ordering time, in seconds
matrix symamd symmmd amdbar mc47bd
pwt 0.78 3.53 0.45 0.44
msc00726
cfd2 7.36 48.68 3.90 3.71
3dtube 6.09 77.75 0.95 0.92
gearbox 16.58 454.49 2.94 2.89
finan512 1.60 2.61 0.92 0.95
Median time
relative to - 6.13 0.39 0.39
symamd

Table

9: Ordering quality
Nonzeros in L (10 3
op. count (10 6 )
matrix symamd symmmd amdbar mc47bd symamd symmmd amdbar mc47bd
pwt 1497 1425 1592 1556 156 140 173 162
msc00726 103 108 111 111 20 22 23 23
cfd2 74832 94444 75008 75008 137470 211328 136476 136476
3dtube 26128 33213 26355 26355 29969 45578 30053 30053
Median result
relative to - 1.26 1.03 1.03 - 1.52 1.08 1.04
colamd
average. It nearly always produced signicantly better orderings. In con-
trast, symamd was always slower than amdbar and mc47bd, although it found
orderings of similar quality (the finan512 is a notable exception, but none
of these methods nds as good an ordering as a tree dissection method [4]).
An ordering algorithm designed for symmetric matrices (amdbar or mc47bd)
is thus superior to one based on a column ordering of M. There may be at
least one important exception, however. In some applications, a better matrix
M is available. Consider an n-by-n nite element matrix constructed as
the summation of e nite elements, where normally e < n. Each nite element
matrix is a dense symmetric submatrix. Suppose element k is nonzero
for all entries a ij for which both i and j are in the set E k . We can construct
a e-by-n matrix M, where row k has the pattern E k . Since M T M has the
same pattern as A, we can compute a column ordering of M and use it for
the Cholesky factorization of A. Colamd would be faster than when using an
M constructed without the knowledge of the nite element structure. The
space to perform the ordering and symbolic analysis is less as well, since M
has fewer nonzeros than A. Although many of the matrices in our test set
arise from nite element problems, only a few small ones are available in
unassembled form, as a collection of nite elements. Thus, we are not able
to evaluate this strategy on large, realistic test matrices.
6

Summary

Two new ordering routines, colamd and symamd, have been presented. For
square nonsymmetric matrices, colamd is much faster and provides better orderings
than Matlab's colmmd routine. It is also faster than symmetric-based
ordering methods (such as amdbar), and uses less storage. For rectangular
matrices (such as those arising in least squares problems and interior point
methods for linear programming), colamd is faster than colmmd and amdbar
and nds orderings of comparable quality. We presented a symmetric ordering
method symamd based on colamd; although it produces orderings as good
as a truly symmetric ordering algorithm (amdbar), it is slower than amdbar.
The colamd and symamd routines are written in ANSI/ISO C, with Matlab-
callable interfaces. Version 2.0 of the code is freely available from the following
sources:
1. University of Florida, http://www.cise.ufl.edu/research/sparse.
2. Netlib, http://www.netlib.org/linalg/colamd/.
3. The MathWorks, Inc., for user-contributed contributions to Matlab,
http://www.mathworks.com. Colamd and symamd are built-in functions
in Matlab Version 6.0.
4. The collected algorithms of the ACM, as Algorithm 8xx, described in
[10].



--R


Compressed graphs and the minimum degree algorithm.
Test collection (non-Hermitian eigenvalue problems)
Solving multistage stochastic programs using tree dissection.


Univ. of Florida sparse matrix collection.


Algorithm 8xx: Colamd
A supernodal approach to sparse partial pivoting.

Distribution of mathematical software via electronic mail.




Yale sparse matrix package
Algorithms and data structures for sparse symmetric Gaussian elimination.

A fast implementation of the minimum degree algorithm using quotient graphs.
An optimal algorithm for symbolic factorization of symmetric matrices.
Computer Solution of Large Sparse Positive De
On the application of the minimum degree algorithm to
An implementation of Gaussian elimination with partial pivoting for sparse systems.
Symbolic factorization for sparse Gaussian elimination with partial pivoting.
Sparse matrices in MATLAB: design and implementation.
Predicting structure in nonsymmetric sparse matrix factorizations.

Sparse partial pivoting in time proportional to arithmetic operations.
Two fast algorithms for sparse matrices: Multiplication and permuted transposition.

Finding good column orderings for sparse QR factorization.
A new polynomial time algorithm for linear program- ming
Approximate de
An approximate minimum degree column ordering algo- rithm

The elimination form of the inverse and its application to linear programming.
An exponentially
Performance of greedy ordering heuristics for sparse Cholesky factorization.
Resende Ramakrishnan and Drezner.
Node selection strategies for bottom-up sparse matrix orderings
On the e

Computing the minimum
Sparse matrix methods for chemical process separation calculations on supercomputers.
Multifrontal vs. frontal techniques for chemical process simulation on supercomput- ers
--TR
A new polynomial-time algorithm for linear programming
Predicting fill for sparse orthogonal factorization
Distribution of mathematical software via electronic mail
Symbolic factorization for sparse Gaussian elimination with partial pivoting
Sparse matrix test problems
A set of level 3 basic linear algebra subprograms
An empirical evaluation of the KORBX algorithms for military airlift applications
A generalized envelope method for sparse factorization by rows
Sparse matrices in matlab
Sparse matrix methods for chemical process separation calculations on supercomputers
Sparsity analysis of the <italic>QR</italic> factorization
An Efficient Algorithm to Compute Row and Column Counts for Sparse Cholesky Factorization
Modification of the minimum-degree algorithm by multiple elimination
Compressed graphs and the minimum degree algorithm
An Approximate Minimum Degree Ordering Algorithm
An Unsymmetric-Pattern Multifrontal Method for Sparse LU Factorization
A combined unifrontal/multifrontal method for unsymmetric sparse matrices
Node Selection Strategies for Bottom-Up Sparse Matrix Ordering
A Supernodal Approach to Sparse Partial Pivoting
Performance of Greedy Ordering Heuristics for Sparse Cholesky Factorization
Two Fast Algorithms for Sparse Matrices: Multiplication and Permuted Transposition
On Algorithms for Obtaining a Maximum Transversal
The Multifrontal Solution of Indefinite Sparse Symmetric Linear
Computer Solution of Large Sparse Positive Definite
Algorithm 836
Algorithm 837

--CTR
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
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
Marco Morandini , Paolo Mantegazza, Using dense storage to solve small sparse linear systems, ACM Transactions on Mathematical Software (TOMS), v.33 n.1, p.5-es, March 2007
Kai Shen, Parallel sparse LU factorization on different message passing platforms, Journal of Parallel and Distributed Computing, v.66 n.11, p.1387-1403, November 2006
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
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
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
Doron Chen , Daniel Cohen-Or , Olga Sorkine , Sivan Toledo, Algebraic analysis of high-pass quantization, ACM Transactions on Graphics (TOG), v.24 n.4, p.1259-1282, October 2005
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
sparse nonsymmetric matrices;ordering methods;linear equations