A combined unifrontal/multifrontal method for unsymmetric sparse matrices.
We discuss the organization of frontal matrices in multifrontal methods for the solution of large sparse sets of unsymmetric linear equations. In the multifrontal method, work on a frontal matrix can be suspended, the frontal matrix can be stored for later reuse, and a new frontal matrix can be generated. There are thus several frontal matrices stored during the factorization, and one or more of these are assembled (summed) when creating a new frontal matrix. Although this means that arbitrary sparsity patterns can be handled efficiently, extra work is required to sum the frontal matrices together and can be costly because indirect addressing is requred. The (uni)frontal method avoids this extra work by factorizing the matrix with a single frontal matrix. Rows and columns are added  to the frontal matrix, and pivot rows and columns are removed. Data movement is simpler, but higher fill-in can result if the matrix cannot be permuted into a variable-band form with small profile. We consider a combined unifrontal/multifrontal algorithm to enable general fill-in reduction orderings to be applied without the data movement of previous multifrontal approaches.  We discuss this technique in the context of a code designed for the solution of sparse systems with unsymmetric pattern.
Introduction
. We consider the direct solution of sets of linear equations
the coefficient matrix A is sparse, unsymmetric, and has a general
nonzero pattern. In a direct approach, a permutation of the matrix A is factorized
into its LU factors, are permutation matrices chosen
to preserve sparsity and maintain numerical accuracy. Many recent algorithms and
software for direct solution of sparse systems are based on a multifrontal approach
[17, 2, 10, 27]. In this paper, we will examine a frontal matrix strategy to be used
within a multifrontal approach.
Frontal and multifrontal methods compute the LU factors of A by using data
structures that permit regular access of memory and the use of dense matrix kernels in
their innermost loops. On supercomputers and high-performance workstations, this
can lead to a significant increase in performance over methods that have irregular
memory access and which do not use dense matrix kernels.
We discuss frontal methods in Section 2 and multifrontal techniques in Section 3,
before discussing the combination of the two methods in Section 4. We consider the
influence of the principle parameter in our approach in Section 5, and the performance
Computer and Information Science and Engineering Department, University of Florida,
Gainesville, Florida, USA. (904) 392-1481, email: davis@cis.ufl.edu. Technical reports and matrices
are available via the World Wide Web at http://www.cis.ufl.edu/~davis, or by anonymous ftp at
ftp.cis.ufl.edu:cis/tech-reports.
y 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. Technical reports, information on HSL, and matrices are
available via the World Wide Web at http://www.cis.rl.ac.uk/struct/ARCD/NUM.html, or by
anonymous ftp at seamus.cc.rl.ac.uk/pub.
T. A. DAVIS AND I. S. DUFF
l
l
l
l
l l
l
Fig. 2.1. Frontal method example
of our new approach in Section 6, before a few concluding remarks and information
on the availability of our codes in Section 7.
2. Frontal methods. In a frontal scheme [15, 26, 31, 32], the factorization
proceeds as a sequence of partial factorizations and eliminations on full submatrices,
called frontal matrices. Although frontal schemes were originally designed for the
solution of finite element problems [26], they can be used on assembled systems [15]
and it is this version that we study in this paper. For general systems, the frontal
matrices can be written as
where all rows are fully summed (that is there is no further contributions to come
to the rows in (2.1)) and the first two block columns are fully summed. This means
that pivots can be chosen from anywhere in the fully summed block comprising the
first two block columns and, within these columns, numerical pivoting with arbitrary
row interchanges can be accommodated since all rows in the frontal matrix are fully-
summed. We assume, without loss of generality, that the pivots that have been chosen
are in the square matrix F 11 of (2.1). F 11 is factorized, multipliers are stored over
F 21 and the Schur complement
using
full matrix kernels. The submatrix consisting of the rows and columns of the frontal
matrix from which pivots have not yet been selected is called the contribution block.
At the next stage, further entries from the original matrix are assembled with the
Schur complement to form another frontal matrix. The overhead is low since each
equation is only assembled once and there is never any assembly of two (or more)
frontal matrices.
An example is shown in Figure 2.1, where two pivot steps have already been
performed on a 5-by-7 frontal matrix (computing the first two rows and columns of
U and L, respectively). Nonzero entries in L and U are shown in lower case. Row 6
has just been assembled into the current 4-by-7 frontal matrix (shown as a solid box).
Columns 3 and 4 are now fully-summed and can be eliminated. After this step, rows 7
and 8 must both be assembled before columns 5 and 6 can be eliminated (the dashed
box, a 4-by-6 frontal matrix containing rows 5 through 8 and columns 5, 6, 7, 8, 9, and
12 of the active submatrix). The frontal matrices are, of course, stored with columns
A COMBINED UNIFRONTAL/MULTIFRONTAL METHOD 3
packed together so no zero columns are held in the frontal matrix delineated by the
dashed box. The dotted box shows the state of the frontal matrix when the next five
pivots can be eliminated. To factorize the 12-by-12 sparse matrix in Figure 2.1, a
working array of size 5-by-7 is required to hold the frontal matrix. Note that,
in

Figure

2.1, the columns are in pivotal order.
One important advantage of the method is that only the frontal matrix need
reside in memory. Entries in A can be read sequentially from disk into the frontal
matrix, one row at a time. Entries in L and U can be written sequentially to disk in
the order they are computed.
The method works well for matrices with small profile, where the profile of a
matrix is a measure of how close the nonzero entries are to the diagonal and is given
by the expression:
a ji 6=0
where it is assumed the diagonal is nonzero so all terms in the summation are non-
negative. If numerical pivoting is not required, fill-in does not increase the profile
U has the same profile as A). To reduce the profile, the frontal method is
typically preceded by an ordering method such as reverse Cuthill-McKee (RCM) [5,
7, 29], which is typically faster than the sparsity-preserving orderings required by a
more general technique like a multifrontal method (such as nested dissection [24] and
minimum degree [1, 25]). A degree update phase, which is typically the most costly
part of a minimum degree algorithm, is not required in the RCM algorithm. However,
for matrices with large profile, the frontal matrix can be large, and an unacceptable
level of fill-in can occur.
3. Multifrontal methods. In a multifrontal scheme for a symmetric matrix
[2, 8, 9, 10, 17, 18, 28], it is normal to use an ordering such as minimum degree to
reduce the fill-in. Such an ordering tends to reduce fill-in much more than profile
reduction orderings. The ordering is combined with a symbolic analysis to generate
an assembly or computational tree, where each node represents the assembly and
elimination operations on a frontal matrix and each edge the transfer of data from child
to parent node. When using the tree to drive the numerical factorization, assembly
and eliminations at any node can proceed as soon as those at the child nodes have been
completed, giving added flexibility for issues such as exploitation of parallelism. As in
the frontal scheme, the complete frontal matrix (2.1) cannot normally be factorized
but only a few steps of Gaussian elimination are possible, after which the remaining
reduced matrix (the Schur complement
needs to be
summed (assembled) with other data at the parent node.
In the unsymmetric-pattern multifrontal method, the tree is replaced by a directed
acyclic graph (dag) [9], and a contribution block may be assembled into more than
one subsequent frontal matrix.
4. Combining the two methods. Let us now consider an approach that
combines some of the best features of the two methods. Assume we have chosen
a pivot and determined a frontal matrix as in a normal multifrontal method. At this
stage, a normal multifrontal method will select as many pivots as it can from the fully
summed part of the frontal matrix, perform the eliminations corresponding to these
pivots, store the contributions to the factors, and store the remaining frontal matrix
for later assembly at the parent node of the assembly tree. The strategy we use in
4 T. A. DAVIS AND I. S. DUFF
our combined method is, after forming the frontal matrix, to hold it as a submatrix
of a larger working array and allow further pivots to be chosen from anywhere within
the frontal matrix. If such a potential pivot lies in the non-fully summed parts of
the frontal matrix then it is necessary to sum its row and column before it could be
used as a pivot. This is possible so long as its fully summed row and column can be
accommodated within the larger working array. In this way, we avoid some of the
data movement and assemblies of the multifrontal method. Although the motivation
is different, the idea of continuing with a frontal matrix for some steps before moving
to another frontal matrix is similar to recent work in implementing frontal schemes
within a domain decomposition environment, for example [23], where several fronts are
used within a unifrontal context. However, in the case of [23], the ordering is done a
priori and no attempt is made to use a minimum degree ordering. Another case where
one continues with a frontal matrix for several assemblies is when relaxed supernode
amalgamation is used [17]. However, this technique delays the selection of pivots and
normally causes more fill-in and operations than our technique. Indeed, the use of our
combined unifrontal/multifrontal approach does not preclude implementing relaxed
supernode amalgamation as well. In comparison with the unifrontal method where
the pivoting is determined entirely from the order of the assemblies, the combined
method requires the selection of pivots to start each new frontal matrix. This implies
either an a priori symbolic ordering or a pivot search during numerical factorization.
The pivot search heuristic requires the degrees of rows and columns, and thus a degree
update phase. The cost of this degree update is clearly a penalty we pay to avoid the
poor fill-in properties of conventional unifrontal orderings.
We now describe how this new frontal matrix strategy is applied in UMFPACK
Version 1.1 [8, 10], in order to obtain a new combined method. However, the new
frontal strategy can be applied to any multifrontal algorithm. UMFPACK does not
use actual degrees but keeps track of upper bounds on the degree of each row and
column. The symmetric analogue of the approximate degree update in UMFPACK
has been incorporated into an approximate minimum degree algorithm (AMD) as
discussed in [1], where the accuracy of our degree bounds is demonstrated.
Our new algorithm consists of several major steps, each of which comprises several
pivot selection and elimination operations. To start a major step, we choose a pivot,
using the degree upper bounds and a numerical threshold test, and then define a
working array which is used for pivoting and elimination operations so long as these
can be performed within it. When this is no longer possible, the working array is
stored and another major step is commenced.
To start a major step, the new algorithm selects a few (by default columns
from those of minimum upper bound degree and computes their patterns and true
degrees. A pivot row is selected on the basis of the upper bound on the row degree
from those rows with nonzero entries in the selected columns. Suppose the pivot row
and column degrees are r and c, respectively. A k-by-l working array is allocated,
typically has a value between 2 and 3, and is fixed
for the entire factorization). The active frontal matrix is c-by-r but is stored in a
k-by-l working array. The pivot row and column are fully assembled into the working
array and define a submatrix of it as the active frontal matrix. The approximate
degree update phase computes the bounds on the degrees of all the rows and columns
in this active frontal matrix and assembles previous contribution blocks into the active
frontal matrix. A row i in a previous contribution block is assembled into the active
frontal matrix if
A COMBINED UNIFRONTAL/MULTIFRONTAL METHOD 5
1. the row index i is in the nonzero pattern of the current pivot column, and
2. if the column indices of the remaining entries in the row are all present in the
nonzero pattern of the current pivot row.
Columns of previous contribution blocks are assembled in an analogous manner.
The major step then continues with a sequence of minor steps at each of which
another pivot is sought. These minor steps are repeated until the factorization can
no longer continue within the current working array, at which point a new major
step is started. When a pivot is chosen in a minor step, its rows and columns are
fully assembled into the working array and redefine the active frontal matrix. The
approximate degree update phase computes the bounds on all rows and columns in
this active frontal matrix and assembles previous contribution blocks into the active
frontal matrix. The corresponding row and column of U and L are computed, but
the updates from this pivot are not necessarily performed immediately. For efficient
use of the Level 3 BLAS it is better to accumulate a few updates (typically 16 or 32,
if possible) and perform them at the same time. To find a pivot in this minor step,
a single candidate column from the contribution block is first selected, choosing one
with least value for the upper bound of the column degree, and any pending updates
are applied to this column. The column is assembled into a separate work vector,
and a pivot row is selected on the basis of upper bound on the row degrees and a
numerical threshold test. Suppose the candidate pivot row and column have degrees
are r 0 and c 0 , respectively. Three conditions apply:
1. If r 0 ? l or c 0 ? k, then factorization can no longer continue within the active
frontal matrix. Any pending updates are applied. The LU factors are stored.
The active contribution block is saved for later assembly into a subsequent
frontal matrix. The major step is now complete.
2. If r then the candidate pivot can fit into the
active frontal matrix without removing the p pivots already stored there.
Factorization continues within the active frontal matrix by
commencing another minor step.
3. Otherwise, if l then the candidate pivot can fit,
but only if some of the previous p pivots are shifted out of the current frontal
matrix. Any pending updates are applied. The LU factors corresponding
to the pivot rows and columns are removed from the front and stored. The
active contribution block is left in place. Set p / 1. Factorization continues
within the active frontal matrix by commencing another minor step.
We know of no previous multifrontal method that considers case 3 (in particular,
UMFPACK Version 1.1 does not). Case 1 does not occur in frontal methods, which are
given a working array large enough to hold the largest frontal matrix (unless numerical
pivoting causes the frontal matrix to grow unexpectedly). Taking simultaneous
advantage of both cases 1 and 3 can significantly reduce the memory requirements and
number of assembly operations, while still allowing the use of orderings that reduce
fill-in.

Figure

4.1 illustrates how the working array is used in UMFPACK Version 1.1
and the combined method. The matrices L 1 , L 2 , U 1 , and U 2 in the figure are the
columns and rows of the LU factors corresponding to the pivots eliminated within
this frontal matrix. The matrix D is the contribution block. The matrix -
U 2 is U 2
with rows in reverse order, and -
columns in reverse order. The arrows
denote how these matrices grow as new pivots are added. When pivots are removed
from the working array in Figure 4.1(b), for case 4 above, the contribution block does
6 T. A. DAVIS AND I. S. DUFF
(a) UMFPACK working array
(b) Combined method working array
empty
empty
U
U
LD
U 121Fig. 4.1. Data structures for the active frontal matrix

Table
Test matrices.
name n jAj sym. discipline comments
gre 1107 1107 5664 0.000 discrete simul. computer system
gemat11 4929 33185 0.001 electric power linear programming basis
migration
lns 3937 3937 25407 0.850 fluid flow linearized Navier-Stokes
shyy161 76480 329762 0.726 fluid flow viscous fully-coupled Navier-Stokes
hydr1 5308 23752 0.004 chemical eng. dynamic simulation
rdist1 4134 94408 0.059 chemical eng. reactive distillation
lhr04 4101 82682 0.015 chemical eng. light hydrocarbon recovery
chemical eng. light hydrocarbon recovery
not need to be repositioned. A similar data organization is employed by MA42 [22].
5. Numerical experiments. We discuss some experiments on the selection of
a value for G in this section and compare the performance of our code with other
sparse matrix codes in the next section.
We show in Table 5.1 the test matrices we will use in our experiments in this
section and the next. The table lists the name, order, number of entries (jAj),
symmetry, the discipline from which the matrix comes, and additional comments.
The symmetry, or more correctly the structural symmetry, is the number of matched
off-diagonal entries over the total number of off-diagonal entries. An entry, a ij (j 6= i),
is matched if a ji is also an entry.
We have performed our experiments on the effect of the value of G on a SUN
SPARCstation 20 Model 41, using the Fortran compiler (f77 version 1.4, with -O4
and -libmil options), and the BLAS from [11] and show the results in Table 5.2. The
table lists the matrix name, the growth factor G, the number of frontal matrices
(major steps), the numerical factorization time in seconds, the total time in seconds,
the number of nonzeros in the LU factors, the memory used, and the floating-point
operation count.
The number of frontal matrices is highest for decreases as G increases
although, because the effect is local, this decrease may not be monotonic. Although
the fill-in and operation count is typically the lowest when the minimum amount of
memory is allocated for each frontal matrix 1), the factorization time is often
high because of the additional data movement required to assemble the contribution
blocks and the fact that the dense matrix kernels are more efficient for larger frontal
matrices.
These results show that our strategy of allocating more space than necessary for
the frontal matrix and choosing pivots that are not from the initial fully summed
A COMBINED UNIFRONTAL/MULTIFRONTAL METHOD 7

Table
Effect of G on the performance of the combined method
Matrix G fronts factor total jL
memory op count
gre 1107 1.0 926 0.51 0.89 0.048 0.160 2.6
2.0 406 0.43 0.74 0.079 0.196 6.7
2.5 323 0.70 1.01 0.100 0.224 9.0
7.0 122 0.75 0.85 0.121 0.235 9.5
2.0 988 0.37 0.70 0.067 0.287 0.8
2.5 873 0.39 0.69 0.077 0.299 1.0
3.0 776 0.33 0.66 0.087 0.312 1.3
5.0 586 0.41 0.70 0.119 0.353 2.0
7.0 538 0.41 0.77 0.148 0.381 2.8
2.0 2324 1.41 3.79 0.122 0.959 7.2
2.5 2326 1.37 3.96 0.122 0.883 7.2
3.0 2322 1.41 3.97 0.123 0.886 8.1
7.0 2280 2.82 5.87 0.153 0.852 16.0
lns 3937 1.0 3208 12.82 17.16 0.474 1.161 95.5
2.0 2011 5.68 9.02 0.494 1.170 84.1
2.5
3.0 1717 7.01 9.29 0.551 1.425 90.0
5.0 1446 6.26 9.30 0.591 1.750 96.4
7.0 1454 11.76 15.14 0.697 2.005 134.1
2.0 3947 0.52 1.57 0.112 0.384 2.7
2.5 3726 0.67 1.46 0.122 0.398 3.3
2.0 306 1.12 1.85 0.278 0.668 10.6
2.5 316 1.21 1.85 0.304 0.694 11.5
3.0 295 1.23 1.69 0.333 0.702 10.4
5.0 160 2.15 2.21 0.493 0.890 14.6
7.0 160 3.17 2.04 0.620 1.038 14.7
2.0 2223 1.90 4.03 0.313 0.824 16.4
2.5 2185 3.09 5.12 0.457 1.051 32.4
5.0 2153 2.51 4.10 0.444 0.995 26.4
7.0 2180 3.70 4.55 0.497 1.195 28.9
block can give substantial gains in execution time with sometimes a small penalty in
storage, although occasionally the storage can be less as well.
6. Performance. In this section, we compare the performance of the combined
unifrontal/multifrontal method (MA38) with the unsymmetric-pattern multifrontal
8 T. A. DAVIS AND I. S. DUFF
method (UMFPACK Version 1.1 [8, 10]), a general sparse matrix factorization
algorithm that is not based on frontal matrices (MA48, [19, 20]), the frontal method
(MA42, [15, 22]), and the symmetric-pattern multifrontal method (MUPS [2]). All
methods can factorize general unsymmetric matrices and all use dense matrix kernels
to some extent [12]. We tested each method on a single processor of a CRAY C-98,
although MUPS is a parallel code. Version 6.0.4.1 of the Fortran compiler (CF77)
was used. Each method (except MA42, which we discuss later) was given 95 Mw of
memory to factorize the matrices listed in Table 5.1. Each method has a set of input
parameters that control its behavior. We used the recommended defaults for most of
these, with a few exceptions that we now indicate.
By default, three of the five methods (MA38, UMFPACK V1.1, and MA48)
preorder a matrix to block triangular form (always preceded by finding a maximum
transversal [14]), and then factorize each block on the diagonal [16]. This can reduce
the work for unsymmetric matrices. We did not perform the preordering, since MA42
and MUPS do not provide these options.
One matrix (lhr71) was so ill-conditioned that it required scaling prior to its
factorization. The scale factors were computed by the Harwell Subroutine Library
routine MC19A [6]. Each row was then subsequently divided by the maximum
absolute value in the row (or column, depending on how the method implements
threshold partial pivoting). No scaling was performed on the other matrices.
By default, MUPS preorders each matrix to maximize the modulus of the smallest
entry on the diagonal (using a maximum transversal algorithm). This is followed by
a minimum degree ordering on the nonzero pattern of A +A T .
For MA42, we first preordered the matrix to reduce its profile using a version
of Sloan's algorithm [30, 21]. MA42 is able to operate both in-core and out-of-
core, using direct access files. It has a finite-element entry option, for which it was
primarily designed. We used the equation entry mode on these matrices, although
some matrices are obtained from finite-element calculations. We tested MA42 both
in-core and out-of-core. The CPU time for the out-of-core factorization was only
slightly higher than the in-core factorization, but the memory usage was much less.
We thus report only the out-of-core results. Note that the CPU time does not include
any additional system or I/O device time required to read and write the direct access
files. The only reliable way to measure this is on a dedicated system. On a CPU-bound
multiuser system, the time spent waiting for I/O would have little effect on
total system throughput.
The symbolic analysis phase in MA42 determines a minimum front size (k-by-l,
which assumes no increase due to numerical pivoting), and minimum sizes of other
integer and real buffers (of total size b, say). We gave the numerical factorization
phase a working array of size 2k-by-2l, and a buffer size of 2b, to allow for numerical
pivoting. In our tables, we show how much space was actually used. The disk space
taken by MA42 is not included in the memory usage.
The results are shown in Table 6.1. For each matrix, the tables list the numerical
factorization time, total factorization time, number of nonzeros in L+U (in millions),
amount of memory used (in millions of words), and floating-point operation count (in
millions of operations) for each method. The total time includes preordering, symbolic
analysis and factorization, and numerical factorization. The time to compute the scale
factors for the lhr71 matrix is not included, since we used the same scaling algorithm
for all methods. For each matrix, the lowest time, memory usage, or operation count
is shown in bold. We compared the solution vectors, x, for each method. We found
A COMBINED UNIFRONTAL/MULTIFRONTAL METHOD 9
that all five methods compute the solutions with comparable accuracy, in terms of
the norm of the residual. We do not give the residual in Table 6.1.
All five codes have factorize-only options that are often much faster than the
combined analysis+factorization phase(s), and indeed the design criterion for some
codes (for example MA48) was to minimize factorize time even if this caused an
increase in the initial analyse time. For MA42, however, the analysis is particularly
simple so normally its overhead is much smaller than for the other codes.
MUPS is shown as failing twice (on shyy161 and lhr71 which are both nearly
singular). In both cases, numerical pivoting caused an increase in storage requirements
above the 95 Mw allocated. MA42 failed to factorize the shyy161 matrix because of
numerical problems. It ran out of disk space during the numerical factorization of
the lhr71 matrix (about 5.8 Gbytes were required). For this matrix, we report the
estimates of the number of entries in the LU factors and the memory usage, as reported
by the symbolic analysis phase of MA42.
Since the codes being compared all offer quite different capabilities and are
designed for different environments, the results should not be interpreted as a
direct comparison between them. However, what we would like to highlight is the
improvements that our new technique brings to the UMFPACK code and that the new
code is at least comparable in performance with other sparse matrix codes. The peak
performance of MA38 is 607 Mflops for the numerical factorization of the psmigr 1
matrix (compared with the theoretical peak of about 1 Gflop for one C-90 processor).
Our new code requires less storage than the original UMFPACK code and sometimes
it requires less than half the storage. In execution time for both analyse/factor and
only, it is generally faster (sometimes by nearly a factor of two) and when
it is slower than the original UMFPACK code it is so by at most about 13% (for the
hydr1 matrix). Over all the codes, MA38 has the fastest analyse+factorize time for
five out of the ten matrices. Except for one matrix (lns 3937) it never takes more than
twice the time of the fastest method. Its memory usage is comparable to the other
in-core methods. MA42 can typically factorize these matrices with the least amount
of core memory. MA42 was originally designed for finite-element entry and so is not
optimized for equation entry. The experiments are run on only one computer which
strongly favors direct addressing and may thus disadvantage MA48 which does not
make such heavy use of the Level 3 BLAS.
7.

Summary

. We have demonstrated how the advantages of the frontal and
multifrontal methods can be combined. The resulting algorithm performs well for
matrices from a wide range of disciplines. We have shown that the combined
unifrontal/multifrontal method gains over our earlier code because it avoids more
indirect addressing and generally performs eliminations on larger frontal matrices.
Other differences between UMFPACK Version 1.1 and the new code (MA38, or
UMFPACK Version 2.0) include an option of overwriting the matrix A with its LU
factors, printing of input and output parameters, a removal of the extra copy of the
numerical values of A, iterative refinement with sparse backward error analysis [4],
more use of Level 3 BLAS within the numerical factorization routine, and a simpler
calling interface. These features improve the robustness of the code and result a
modest decrease in memory usage.
The combined unifrontal/multifrontal method is available as the Fortran 77 codes,
UMFPACK Version 2.0 in Netlib [13], 1 and MA38 in Release 12 of the Harwell
1 UMFPACK Version 2.0, in Netlib, may only be used for research, education, or benchmarking
T. A. DAVIS AND I. S. DUFF
Subroutine Library [3]. 2
8.

Acknowledgements

. We would like to thank Nick Gould, John Reid, and
Jennifer Scott from the Rutherford Appleton Laboratory for their helpful comments
on a draft of this report.



--R

An approximate minimum degree ordering algorithm
Vectorization of a multiprocessor multifrontal code

Solving sparse linear systems with sparse backward error
A linear time implementation of the reverse Cuthill-Mckee algorithm
On the automatic scaling of matrices for Gaussian elimination
Reducing the bandwidth of sparse symmetric matrices
Users' guide to the unsymmetric-pattern multifrontal package (UMFPACK




Distribution of mathematical software via electronic mail
On algorithms for obtaining a maximum transversal

An implementation of Tarjan's algorithm for the block triangularization of a matrix



The design of MA48
The use of profile reduction algorithms with a frontal code

The use of multiple fronts in Gaussian elimination
Computer Solution of Large Sparse Positive Definite Systems

A frontal solution program for finite element analysis
The multifrontal method for sparse matrix solution: Theory and Practice
The multifrontal method for sparse matrix solution: Theory and practice
Comparative analysis of the Cuthill-Mckee and the reverse Cuthill-Mckee ordering algorithms for sparse matrices
An algorithm for profile and wavefront reduction of sparse matrices
frontal techniques for chemical process simulation on supercomputers
Supercomputing strategies for the design and analysis of complex separation systems
--TR
Distribution of mathematical software via electronic mail
Sparse matrix test problems
The evolution of the minimum degree ordering algorithm
A set of level 3 basic linear algebra subprograms
The multifrontal method for sparse matrix solution
Sparse matrix methods for chemical process separation calculations on supercomputers
The design of a new frontal code for solving sparse, unsymmetric systems
The design of MA48
Stable finite elements for problems with wild coefficients
An Approximate Minimum Degree Ordering Algorithm
An Unsymmetric-Pattern Multifrontal Method for Sparse LU Factorization
The Multifrontal Solution of Indefinite Sparse Symmetric Linear
Computer Solution of Large Sparse Positive Definite
Reducing the bandwidth of sparse symmetric matrices
A Supernodal Approach to Sparse Partial Pivoting

--CTR
Eero Vainikko , Ivan G. Graham, A parallel solver for PDE systems and application to the incompressible Navier-Stokes equations, Applied Numerical Mathematics, v.49 n.1, p.97-116, April 2004
J.-R. de Dreuzy , J. Erhel, Efficient algorithms for the determination of the connected fracture network and the solution to the steady-state flow equation in fracture networks, Computers & Geosciences, v.29 n.1, p.107-111, February
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
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
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
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
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
J. Candy , R. E. Waltz, An Eulerian gyrokinetic-Maxwell solver, Journal of Computational Physics, v.186
A. Guardone , L. Vigevano, Finite element/volume solution to axisymmetric conservation laws, Journal of Computational Physics, v.224 n.2, p.489-518, June, 2007
frontal methods;sparse unsymmetric matrices;linear equations;multifrontal methods