An Unsymmetrized Multifrontal LU Factorization.
A well-known approach to computing the LU factorization of a general unsymmetric matrix A is to build the elimination tree associated with the pattern of the symmetric matrix A + AT and use it as a computational graph to drive the numerical factorization. This approach, although very efficient on a large range of unsymmetric matrices, does not capture the unsymmetric structure of the matrices. We introduce a new algorithm which detects and exploits the structural asymmetry of the submatrices involved during the processing of the elimination tree. We show that with the new algorithm, significant gains, both in memory and in time, to perform the factorization can be obtained.
Introduction
We consider the direct solution of sparse linear equations based on a multifrontal approach.
The systems are of the form A is an n  n unsymmetric sparse matrix. The
multifrontal method has been developed by Du and Reid [11, 12] for computing the solution
of indenite sparse symmetric linear equations using Gaussian elimination and then has been
extended to solve more general unsymmetric matrices by Du and Reid [13].
The multifrontal method belongs to the class of methods that separate the factorization into
an analysis phase and a numerical factorization. The analysis phase involves a reordering
step, which will reduce the ll-in during numerical factorization and a symbolic phase that
This work was supported by the Director, Oce of Science, Division of Mathematical, Information, and
Computational Sciences of the U.S. Department of Energy under contract number DE-AC03-76SF00098.
y amestoy@enseeiht.fr, ENSEEIHT-IRIT, 2 rue Camichel 31071 Toulouse (France) and PRamestoy@lbl.gov,
NERSC, Lawrence Berkeley National Laboratory, 1 Cyclotron Rd. Berkeley CA 94720
z NERSC, Lawrence Berkeley National Laboratory, 1 Cyclotron Rd., Berkeley CA 94720
builds the computational tree, so called elimination tree [9, 18, 20], whose structure gives the
dependency graph of the multifrontal approach. The analysis phase is generally not concerned
with numerical values and is only based on the sparsity pattern of the matrix.
As far as the analysis phase is concerned, the approaches introduced by Du and Reid
for both symmetric and unsymmetric matrices are almost identical. When the matrix is
unsymmetric, the structurally symmetric matrix the summation is
performed symbolically, is used in place of the original matrix A. The elimination tree of
the unsymmetric LU factorization is thus identical to that of the Cholesky factorization of the
symmetrized matrix M.
To control the growth of the factors during LU factorization, partial threshold pivoting is
used during the numerical factorization phase. The pivot order, used during the analysis to
build the elimination tree might not be respected. Numerical pivoting can then result in an
increase in the estimated size of the factors and in the number of operations. To improve the
numerical behaviour of the multifrontal approach it is common to involve a step of preprocessing
based on the numerical values. In fact if the matrix is not well-scaled, which means that the
entries in the original matrix do not have the same order of magnitude, a good prescaling
of the matrix can have a signicant impact on the accuracy and performance of the sparse
solver. In some cases it is also very benecial to precede the ordering by performing an
unsymmetric permutation to place large entries on the diagonal. Du and Koster [10] have
designed algorithms to permute large entries onto the diagonal and have shown that it can very
signicantly improve the behaviour of multifrontal solvers.
The multifrontal approach by Du and Reid [13] is used in the Harwell Subroutine Library code
3] and in the distributed memory code MUMPS developed in the context of the PARASOL
project (EU ESPRIT IV LTR project 20160) [4, 5]. Another way to represent the symbolic LU
factorization of a structurally unsymmetric matrix is to use directed acyclic graphs (see for
example [14, 15]). These structures more costly and complicated to handle than a tree, capture
better the asymmetry of the matrix. Davis and Du [6] implicitly use this structure to drive
their unsymmetric-pattern multifrontal approach.
We explain, in this article, how to use the simple elimination tree structure of the symmetric
matrix M to detect, during the numerical factorization phase, structural asymmetry in the
factors. We show that, with the new factorization phase, we very signicantly reduce the
computational time, the size of the LU factors and the total memory requirement with respect
to the standard multifrontal approach [13]. In Section 2, we rst recall the main properties of
the elimination tree and describe the standard multifrontal factorization algorithm. We then
introduce the new algorithm and use a simple example to show the benets that can be expected
from the new approach. In Section 3, our set of test matrices is introduced. We analyse the
performance gains (in terms of size of the factors, memory requirement and factorization time)
of the new approach with respect to the standard multifrontal code on our set of test matrices.
We add some concluding remarks in Section 4
2 Description of the multifrontal factorization algorithms
Let A be an unsymmetric matrix and let M denote the structurally symmetric matrix A+A T .
The elimination tree is dened using the structure of the Cholesky factors of M. If the matrix
M is reducible then the tree will be a forest. Liu [18] denes the elimination tree as the
transitive reduction of the directed graph of the Cholesky factors of M. The characterization
of the elimination tree and the description of its properties are beyond the scope of this article.
In our context, we are interested in the elimination tree only as the computational graph for
the multifrontal factorization. For a complete description of the elimination tree the reader can
In the multifrontal approaches, we actually use an amalgamated elimination tree, referred to as
the assembly tree [12] which can be obtained from the classical elimination tree. Each node
of the assembly tree corresponds to Gaussian elimination operations on a full submatrix, called
a frontal matrix. The frontal matrix can be partitioned as shown in Figure 1.
fully summed rows -
partly summed rows -
fully summed columns
partly summed columns

Figure

1: Partitioning of a frontal matrix.
Each frontal matrix factorization involves the computation of a block of columns of L,
termed fully summed columns of the frontal matrix, a block of rows of U, termed fully
summed rows, and the computation of a Schur complement matrix F 22 F 21 F 1
called a
contribution block. The rows (columns) of the F 22 block are referred to as partly summed
rows (columns).
original nonzero
new entry in M

Figure

2: Example of matrix A and
The unsymmetric matrix A, on the left-hand side of Figure 2, will be used to illustrate the
main properties of the assembly tree and to introduce the new algorithm. In Figures 2 and 3
an \X" denotes a nonzero position from the original matrix A and a \ i
" corresponds to a
new entry introduced during symmetrization. In Figure 3, we indicate the structure of the lled
matrix is the matrix of the Cholesky factor of M. Entries with an \F"
corresponds to ll-in entries in the L factor.
F
F
F
original nonzero
Fill-in (w.r.t. M)
new entry in M

Figure

3: Structure of the Cholesky factors of the matrix M.
The matrix M F is used to dene the assembly tree (see Figure 4) associated with the multifrontal
LU factorization of the matrix A. From the fact that the factorization is based on the assembly
tree associated with the Cholesky factorization of M F , it results that
where Struct() denotes the matrix pattern. Let us denote by structural zero a numerical
zero that does not result from numerical cancellation. Typically, due to the symmetrization,
the matrix M might contain many structural zeros that will propagate during the numerical
factorization phase. What has motivated our work is the following question. Is it possible,
during the processing of the assembly tree to eciently detect and remove structural zeros that
appear in matrix M F and that are direct or indirect consequence of the symmetrization of
Although it is not so clear from the structure of the matrix M F , we will show that
blocks of structural zeros can be identied during the processing of the assembly tree.
In the following, we rst describe how the assembly tree is exploited during the standard
multifrontal algorithm. We then report and analyse the sparsity structure of the frontal matrices
involved in the processing of the assembly tree associated with our example matrix. Based on
these observations, we will introduce the new factorization algorithm.
The assembly tree is rooted (a node of the tree called the root is chosen to give an orientation
to the tree) and is processed from the leaf nodes to the root node. If two nodes are adjacent
in the tree, then the one nearer the root is the parent node, and the other is termed its
child. Each edge of the assembly tree indicates a data dependency between parent and child.
It involves sending a contribution block from the child to the parent. A parent node process
will start when the processes associated with all of its children are completed.
U
U
U
U
original nonzero
arrowhead of var. 3
contribution block
Root
(2)
(1)

Figure

4: Assembly tree associated with our test matrix.
For example, in Figure 4, node (3) must wait for the completion of nodes (1) and (2) before
starting its computations. The subset of variables which can be used as pivots (boldface
variables in Figure 4) are the fully summed variables of node (k). The contribution blocks
of the children and the entries from the original matrix corresponding to the fully summed
variables of node (k) are used to build the frontal matrix of the node. This will be referred
to as the assembly process. During the assembly process of a frontal matrix, we need for
each fully summed variable j, to access the nonzero elements in the original matrix that are in
rows/columns of indices greater than j. A way to eciently access the original matrix is to store
it in arrowheads according to the reordered matrix. For example during the assembly process
of node (3) the arrowheads of variables 3 and 4 from matrix A together with the contribution
blocks of nodes (1) and (2) are used to assemble the frontal matrix of node (3). One should
note that, by construction, the list of indices in the partly summed rows is identical to that
of the partly summed columns (row and column indices of block F 22 in Figure 1). Therefore,
during the assembly process, only the list of row indices of the partly summed rows is built. This
list is obtained by merging all the row and column indices of the arrowheads of the matrix A
with the row indices of the contribution blocks of all the sons. Once the structure of the frontal
matrix is built, the numerical values from both the arrowheads and the contribution blocks can
be assembled at the right place in the frontal matrix. The
oating point operations involved
during the assembly process will be referred to as assembly operations (only additions)
whereas
oating-point operations involved during the factorization of the frontal matrices will
be referred to as elimination operations.
Partial threshold pivoting is used to control the element growth in the factors. Note that
pivots can be chosen only from within the block F 11 of the frontal matrix. The LU factors
corresponding to the fully summed variables are computed and a new contribution block is
produced. When a fully summed variable of node (k) cannot be eliminated during the node
process because of numerical considerations, then the corresponding arrowhead in the frontal
matrix is added to the contribution block and the fully summed variable will be included in the
fully summed variables at the parent of node (k). This process creates additional ll-in in the
LU factors.
In a multifrontal algorithm, we have to provide space for the frontal matrices and the
contribution blocks, and to reserve space for storing the factors. We need working space to
store both real and integer information. This will be referred to as the total working space
of the factorization phase. The same integer array can be used to describe a frontal matrix, its
corresponding LU factors and its contribution block. The management of the integer working
array can thus be done in a simple and ecient way. In a uniprocessor environment, it is possible
to determine the order in which the assembly tree will be processed. Furthermore, if we process
the assembly tree with a depth rst search order, we can use a stack to manage the storage
of the factors and the contribution blocks. This mechanism is ecient both in terms of total
memory requirement and amount of data movement (see [12]). A stack mechanism, starting
from the beginning of the real working array, is used to store the LU factors. Another stack
mechanism starting from the end of the real working array is used to store the contribution
blocks. After the assembly phase of a node the working space used by the contribution blocks
of its children can be freed and, because the assembly tree is processed with a depth rst search
order, the contribution blocks will always be at the top of the stack. In the remainder of this
paper, the maximum stack size of the contribution blocks will be referred to as the maximum
stack size.
The standard and new algorithms for multifrontal factorization
During a multifrontal factorization, each frontal matrix can be viewed as the minimum structure
to perform the elimination of the fully summed variables and to carry the contribution blocks
from all of its sons. In Figure 5, we have a closer look at the frontal matrices involved in the
processing of the assembly tree of Figure 4 to identify the structural zeros.
We report, beside each node, the structure of the factorized frontal matrix assuming that the
pivots are chosen down the diagonal of the fully summed block and in order (i.e. no numerical
pivoting is required). An \X" corresponds to a nonzero entry and a \O" corresponds to a
structural zero.
One can see that, for our example, the frontal matrices have many structural zeros. There
are two kinds of structural zeros: those forming a complete zero column (or row), and more
isolated zero entries in a nonzero column or row (for example entries (4,3) and (4,7) in the
e
standard algorithm.
frontal matrix of node (3)). If one knows how to detect a partly summed row (or column)
with only structural zeros then the corresponding row (or column) can be suppressed from the
frontal matrix because this row (or column) will not add any contribution to the father node.
Structural zero rows (or columns) can be detected during the assembly process of a frontal
matrix because of the following property: if a row (or column) index does not appear in the
row (or column) indices both of the arrowheads of the original matrix and of the contribution
blocks of the sons, then this index will correspond to a row (or column) with only structural
zeros. This property is used to deduce the assembly process of the new algorithm. Note that
if the matrix is not structurally decient then each fully summed row (or column) must have
at least one nonzero entry. Therefore, we can restrict our search for zero rows (columns) to the
partly summed rows (columns).
In the new assembly algorithm, the list of indices of the partly summed rows of a frontal matrix
is dened as the merge of the row indices in the arrowheads of the fully summed variables of the
node with the row indices of the contribution blocks of its sons. The column indices are dened
similarly. As it is illustrated in Figure 6, the new assembly process can result in signicant
modications in the processing of the assembly tree. For example, on node (1), row 3 and
column 5 are suppressed from the frontal matrix; on node (2), all the partly summed rows are
suppressed; on node (3), row 7 and column 5 are suppressed. As it can be noticed in Figure 6,
frontal matrices naturally become unsymmetric in structure.
We nally indicate in Figure 7 the structure of the LU factors obtained with the new algorithm.

Figure

Processing the assembly tree associated with the matrix A in Figure 2 using the new
algorithm.
This should be compared to the matrix M F in Figure 3 showing the structure of the factors
obtained with the standard algorithm. It can be seen that nonzero entries corresponding to
original nonzero
Fill-in (w.r.t. M)
new entry in M

Figure

7: Structure of the LU factors obtained with the new algorithm.
ll-in (for example (7,4) in M F ) or introduced during the symmetrization of M (for example
in M F ) might be suppressed by the new algorithm. On the other hand, the new algorithm
will never suppress structural zeros in a block of fully summed variables (for example (4,3) in
node (3) of Figure 5). On our small example, the total number of entries in the factors reduces
from 31 to 23.
Comparing Figures 5 and 6, one can notice that the new algorithm might also lead to a signicant
reduction in both the number of operations involved during the assembly process and the
maximum stack size. The latest combined with a reduction in the size of the factors will result
in a reduction in the total working space. On our example the number of assembly operations
drops from 30 to 20 (18 entries from A plus 2 from the contribution blocks). The maximum
stack size reduces from 8 to 1 (obtained in both cases after stacking the contribution blocks of
nodes (1) and (2)).
3 Results and performance analysis
We describe in Table 1 the set of test matrices (order, number of nonzero entries, structural
symmetry and origin). We dene the structural symmetry as the percentage of the number
of nonzeros matched by nonzeros in symmetric locations over the total number of entries. A
symmetric matrix has a value of 100. Although, our performance analysis will focus on matrices
with a relatively small structural symmetry, all classes of unsymmetric matrices are represented
in this set. The selected matrices come from the forthcoming Rutherford-Boeing Sparse Matrix
Collection [8] 1 , Tim Davis collection 2 , and SPARSEKIT2 3 .
The Harwell Subroutine Library [16] code ma41 has been used to obtain the results for the
standard multifrontal method. The factorization phase of ma41 has then been modied with the
new algorithm. The ma41 code has a set of parameters to control its eciency. We have used the
default values for our target computer. Approximate minimum degree ordering [1] has been used
to reorder the matrix. As we have mentioned in the Introduction, it is often quite benecial for
very unsymmetric matrices to precede the ordering by performing an unsymmetric permutation
to place large entries on the diagonal and then scaling the matrix so that the diagonal entries
are all of modulus one and the o-diagonals have modulus less than or equal to one. We use
the Harwell Subroutine Library code mc64 [10] to perform this preordering and scaling on all
matrices of structural symmetry smaller than 55. When mc64 is not used, our matrices are
always row and column scaled (each row/column is divided by its maximum value). All results
presented in this section, have been obtained on one processor (R10000 MIMPS RISC 64-bit
processor) of the SGI Cray Origin 2000 from Parallab (University of Bergen, Norway). The
processor runs at a frequency of 195 Mhertz and has a peak performance of 400 M
ops per
second.
Web page http://www.cse.clrc.ac.uk/Activity/SparseMatrices/
Web page http://www.cise.ufl.edu/ davis/sparse/
3 Web page http://iftp.cs.umn.edu/pub/sparse/
Matrix name Order No. entries StrSym Origin (Discipline)
av4408 4408 95752 0 Vavasis (Partial di. eqn.) [21]
bbmat 38744 1771722 54 Rutherford-Boeing (CFD)
cavity26 4562 138187 95 SPARSEKIT2 (CFD)
ex11 16614 1096948 100 SPARSEKIT2 (CFD)
goodwin 7320 324784 64 Davis (CFD)
lhr14c 14270 307858 1 Davis (Chemical engineering)
lhr17c 17576 381975 0 Davis (Chemical engineering)
lhr34c 35152 764014 0 Davis (Chemical engineering)
lhr71c 70304 1528092 0 Davis (Chemical engineering)
lns 3937 3937 25407 87 Rutherford-Boeing (CFD)
Rutherford-Boeing (Economics)
Rutherford-Boeing (Demography)
raefsky6 3402 137845 2 Davis (Structural engineering)
Rutherford-Boeing (Chemical engineering)
rim 22560 1014951
sherman5 3312 20793 78 Rutherford-Boeing (Oil reservoir simul.)
shyy161 76480 329762 77 Davis (CFD)
twotone 120750 1224224 28 Davis (Circuit simulation)
wang4 26068 177196 100 Rutherford-Boeing (Semiconductor)

Table

1: Test matrices. StrSym denotes the structural symmetry.
In the following graphs, we report the performance ratios of the new factorization algorithm
over the standard algorithm. Matrices are sorted by increasing structural symmetry of the
matrix to be factored, i.e. after application of the column permutation when mc64 is used. We
use the same matrix order in the graphs and in the complete set of results provided in Tables 2
and 3. In this way, one can easily nd, given a point in the graph, its corresponding entry in
the tables.
On the complete set of test matrices, we rst analyse in Figure 8 what is probably of main
concern for the user of a sparse solver, i.e. the time to factor the matrix and the total working
space (as dened in the previous section). In Figure 8, we divide the matrices into three
categories: matrices of structural symmetry smaller than 50 for which the time reduction is
between 20% and 80%, matrices whose structural symmetry is between 50 and 80 for which the
time reduction is between 3% and 20% and nearly structurally symmetric matrices for which
there is almost no dierence between the standard and new version. It is interesting to notice
that even on symmetric matrices the added work to detect asymmetry does not aect much the
performance of the factorization. In the remainder of this paper, we will not report results on
almost structurally symmetric matrices (symmetry greater than 80).
43 43 48
new
over
standard
version
total working space

Figure

8: Study of the factorization time and the total working space. sy (in abscissa)
corresponds to structurally symmetric matrices.
In

Figure

9, we relate the gain in the factorization time with the reduction in the number of
elimination operations and in the number of assembly operations. Although the number of
operations due to the assembly is always much smaller that the number of operations involved
during factorization (see Tables 2 and 3), the assembly process can still represent an important
part of the time spent in the factorization phase (see for example [2]). This is illustrated in

Figure

9 where we see that the high reduction in the number of assembly operations (more
than 50%) signicantly contributes to reducing the factorization time. Note that on a relatively
large matrix (twotone) of symmetry 57 still signicant gains in time and in number of assembly
operations (more than 40%) can be obtained. In Figure 10, we relate the total working space
reduction to the size of the factors and to the maximum stack size. Although a reduction in the
maximum size of the stack might not always introduce a reduction in the total working space,
we see that in practice it is often the case. An extreme example of this reduction is matrix
orani678 of symmetry 9 (see Table 2) for which maximum stack size is reduced by more than
one order of magnitude (5482454 to 457312). Finally, we notice that a large reduction in the
maximum stack size (Figure 10) will generally correspond to a large reduction in the number
of assembly operations (Figure 9).
Structural Symmetry
new
over
standard
version
oper for factors
oper for assem.

Figure

9: Impact of the reduction in the number of
oating point operations on the time.
43 43 48
new
over
standard
version
maximum size of stack
space for LU factors
total working space

Figure

10: Correlation between the factor size, the maximum stack size, and the total working
space.
Matrix Total Operations in Facto.
(Str.Sym.) Version LU stack space Elimin. Assemb. Time
raefsky6 Stnd 1509016 606458 2017106 4.795E+08 4.348E+06 2.17
(2) New 998064 145575 1119135 2.313E+08 1.049E+06 1.14
raefsky5 Stnd 1757680 378792 2095082 3.746E+08 3.874E+06 1.83
av4408 Stnd 551354 227787 677254 6.872E+07 1.451E+06 0.46
lhr14c Stnd 2167304 415090 2439502 2.092E+08 7.944E+06 1.77
lhr34c Stnd 5613656 755710 6249187 6.282E+08 2.064E+07 5.24
lhr17c Stnd 2813418 639172 3204801 3.089E+08 1.085E+07 2.65
lhr71c Stnd 11615170 729857 12711920 1.402E+09 4.304E+07 13.16
twotone Stnd 22086166 15899616 34489449 2.933E+10 2.171E+08 183.84
onetone1 Stnd 4713485 3348215 6212037 2.282E+09 2.675E+07 14.54
rdist1 Stnd 258096 53767 279999 8.150E+06 5.054E+05 0.13

Table

2: Comparison of the standard (Stnd) and the new algorithms on matrices of structural
symmetry < 50.
Matrix Total Operations in Facto.
(Str.Sym.) Version LU stack space Elimin. Assemb. Time
bbmat Stnd 44111480 8351266 48035816 3.676E+10 2.283E+08 185.75
utm3060 Stnd 324640 78679 806970 2.683E+07 6.973E+05 0.22
utm5940 Stnd 701496 131839 2799224 6.640E+07 1.529E+06 0.51
onetone2 Stnd 2253553 898540 385816 5.085E+08 7.628E+06 3.88
goodwin Stnd 1264140 307777 1385604 1.612E+08 2.841E+06 1.05
rim Stnd 4127204 833290 4371615 5.648E+08 9.194E+06 3.37
shyy161 Stnd 7437816 377535 290589 9.945E+08 1.178E+07 6.56
shyy41 Stnd 251336 28523 8137032 1.036E+07 6.337E+05 0.14
sherman5 Stnd 167412 61227 217769 1.284E+07 4.414E+05 0.13
lns 3937 Stnd 285517 89578 335285 1.920E+07 5.482E+05 0.20
cavity15 Stnd 202629 33004 230111 1.033E+07 3.453E+05 0.10
cavity26 Stnd 394164 58589 447293 2.433E+07 6.877E+05 0.22
ex11 Stnd 11981558 3960507 13614400 6.678E+09 3.836E+07 27.76
fidapm11 Stnd 15997220 4863371 19069150 9.599E+09 4.705E+07 39.78
olaf1 Stnd 5880174 1506068 6794919 1.965E+09 1.685E+07 8.91
wang4 Stnd 11561486 5063375 15900237 1.048E+10 4.087E+07 46.19

Table

3: Comparison of the standard (Stnd) and the new algorithms on matrices of structural
symmetry >= 50
Concluding remarks
We have described a modication of the standard multifrontal LUfactorization algorithm that
can lead to a signicant reduction in both the computational time and the total working space.
The standard multifrontal algorithm [13] for unsymmetric matrices is based on the assembly
tree of a symmetrized matrix and involves frontal matrices symmetric in structure. Therefore, it
produces LU factors such that the matrix F = L+U is symmetric in structure. This approach
is currently used in the context of two publically available packages (ma41 [2, 3] and MUMPS [5, 4])
and has the advantage, with respect to other unsymmetric factorization algorithms [6, 7, 17],
of having the LU factorization based on the processing of an assembly tree, while the other
approaches explicitly or implicitly use a more complex to handle graph structure.
We have demonstrated that, based on the same assembly tree, one can derive a new multifrontal
algorithm that will introduce asymmetry in the frontal matrices and in the matrix of the factors
F. The detection of the asymmetry is only based on structural information and is not costly to
compute as it has been illustrated on structurally symmetric matrices, for which both algorithms
behave similarly. On a set of unsymmetric matrices, we have shown that the new algorithm will
reduce both the factor size and the number of operations by a signicant factor. We have also
observed that the reduction in the number of indirect memory access operations involved during
the assembly process is generally much higher than the reduction in the number of elimination
operations. Finally, we have noticed that the gain in the maximum stack size is also relatively
high and is comparable to the gain in the number of assembly operations.
Space for Operations
LU Stack Total Elim Assemb. Time
mean 0.79 0.35 0.69 0.67 0.41 0.59
median 0.80 0.35 0.76 0.68 0.40 0.61
50  Structural symmetry < 80
mean 0.93 0.85 0.93 0.89 0.82 0.86
median 0.95 0.91 0.95 0.93 0.89 0.90

Table

4: Performance ratios of the new algorithm over the standard algorithm.
To conclude, we report in Table 4 a summary of the results (mean and median) obtained on
the test matrices with structural symmetry smaller than 80. For very unsymmetric matrices
(structural symmetry smaller that 50), we obtain an average reduction of 31% in the total
working space and of 41% in the factorization time. The maximum stack size and the number
of assembly operations are reduced by respectively 65% and 59%. Finally, it is interesting to
observe that, even on fairly symmetric matrices (50  structural symmetry < 80), it can still
be worth trying to identify and exploit asymmetry during the processing of the assembly tree.

Acknowledgements

. We want to thank Horst Simon and Esmond Ng who gave us the
opportunity to work at NERSC (LBNL) for one year. We are grateful to Sherry Li for helpful
comments on an early version of this paper.



--R














Exploiting structural symmetry in unsymmetric sparse symbolic factorization.
Elimination structures for unsymmetric sparse lu factors.

A scalable sparse direct solver using static pivoting.
The role of elimination trees in sparse factorization.
The multifrontal method for sparse matrix solution: Theory and Practice.
A new implementation of sparse Gaussian elimination.

--TR

--CTR
Gianmarco Manzini , Mario Putti, Mesh locking effects in the finite volume solution of 2-D anisotropic diffusion equations, Journal of Computational Physics, v.220 n.2, p.751-771, January, 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
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
Wook Ryol Hwang , Martien A. Hulsen , Han E. H. Meijer, Direct simulation of particle suspensions in sliding bi-periodic frames, Journal of Computational Physics, v.194 n.2, p.742-772, March 2004
See Jo Kim , Wook Ryol Hwang, Direct numerical simulations of droplet emulsions in sliding bi-periodic frames using the level-set method, Journal of Computational Physics, v.225 n.1, p.615-634, July, 2007
elimination tree;gaussian elimination;unsymmetric matrices;multifrontal methods;sparse linear equations