Improved Symbolic and Numerical Factorization Algorithms for Unsymmetric Sparse Matrices.
We present algorithms for the symbolic and numerical factorization phases in the direct solution of sparse unsymmetric systems of linear equations. We have modified a classical symbolic factorization algorithm for unsymmetric matrices to inexpensively compute minimal elimination structures.   We give an efficient algorithm to compute a near-minimal data-dependency graph for unsymmetric multifrontal factorization that is valid irrespective of the amount of dynamic pivoting performed during factorization.  Finally, we describe an unsymmetric-pattern multifrontal algorithm for Gaussian elimination with partial pivoting  that uses the task- and data-dependency graphs computed during the symbolic phase. These algorithms have been implemented in WSMP---an industrial strength sparse solver package---and have enabled WSMP to significantly outperform other similar solvers.  We present experimental results to demonstrate the merits of the new algorithms.
Introduction
. Typical
dire ct
solve rs for
ge ne ral
sparse syste ms of
line ar
e quations of
the form
have four distinct
phase s: analysis, comprising orde
ring for fill-in
re duction and symbolic factorization;
nume rical factorization of
the sparse coe #cie nt matrix A into triangular factors L and U using
Gaussiane limina-
tion with partial pivoting; forward and
backwarde limination to
solve for x using
the triangular factors L and U and
the right-hand-side ve ctor b; and
ite rative re fine me nt
of
the compute d solution. In this
pape r,
we de scribe some of
the algorithms that
are use d in
the unsymme tric symbolic and
nume rical factorization
phase s of
the Watson
Sparse Matrix
Package (WSMP)-a
high-pe rformance and robust
software for solving
ge ne ral
sparse line ar
syste ms.
The se algorithms
are crucial to WSMP's
pe rformance ,
which has
be e n shown to
be significantly
be tte r than that of
othe r similar
solve rs [18].
An important contribution of this
pape r is to show that, contrary to
conve ntional
wisdom, it is
possible to symbolically
de te rmine a static communication
patte rn for
paralle l
unsymme tric
sparse LU
ve n in
the pre se nce of partial pivoting.
The proce ss of factoring a
sparse matrix can
be e xpre sse d by a
dire cte d acyclic
task-de pe nde ncy graph (task-DAG).
The ve rtice s of this DAG
corre spond to
the tasks
of factoring rows or columns, or groups of rows and columns, of
the sparse matrix,
and
the e dge s
corre spond to
the de pe nde ncie s
be twe e n
the tasks. A task is
re ady
fore xe cution if and only if all tasks with
incominge dge s to it
have comple te d. In
addition to a task-DAG,
the re is a
data-de pe nde ncy graph (data-DAG)
associate d
with
sparse matrix factorization.
The ve rte x
se t of
the data-DAG is
the same as that
of
the task-DAG for a
give n
sparse matrix.
Ane dge from a
ve rte x i to a
ve rte x j in
the data-DAG
de note s that at
le ast
some of
the output data of task i is
re quire d as
input by task j. In this
pape r,
we de fine task i as
the task of computing column i of
L and row i of U .
Once the tasks
are de fine d,
the task-DAG is
unique to a
sparse # Received by the editors October 3, 2001; accepted for publication (in revised form) by E. G. Ng
October 18, 2002; published electronically December 3, 2002.
http://www.siam.org/journals/simax/24-2/39604.html
T.J.W atson Research Center, P.O. Box 218, Yorktown Heights, NY 10598 (anshul@watson.
ibm.com).
matrix for a
give n
pe rmutation of rows and columns;
howe ve r,
the data-DAG is a
function of
the sparse factorization algorithm. Multifrontal algorithms [9, 14, 23] for
sparse factorization can work with a minimal data-DAG
(i.e ., a data-DAG with
the smalle st
possible
dge s) for a
give n matrix.
In
the case of
symme tric
sparse matrice s,
the minimal task- and data-DAGs for
the factorization
proce ss
are a
tre e calle d
the e limination
tre e [22].
Howe ve r, for un-
symme tric
sparse matrice s,
the task- and data-DAGs
are ge ne ral DAGs.
More ove r,
the e dge -se t of
the minimal data-DAG for
unsymme tric
sparse factorization can
be a
rse t of
the e dge -se t of a task-DAG.
Gilbe rt and Liu [16]
de scribe e limination
structure s for
unsymme tric
sparse LU factors and
give an algorithm for
sparse unsym-
me tric symbolic factorization.
The se e limination
structure s
are two DAGs that
are transitive re ductions of
the graphs of
the factor
matrice s L and U ,
re spe ctive ly, and
can
be use d to
de rive a task-DAG for
sparse LU factorization.
Some re se arche rs
have argue d that computing
ane xact
transitive re duction can
be tooe xpe nsive [9, 15] and
have propose d using subminimal DAGs with
more e dge s than
ne ce ssary.
Howe ve r,
trave rsing
unne ce ssary
dge s during
nume rical factorization can
be a
source of
ove rhe ad.
More ove r, in a
paralle l
imple me ntation,e xtra
dge s can
be pote ntial
source s of
unne ce ssary synchronization or communication.
In this
pape r,
we show how a
re lative ly straightforward modification to
Gilbe rt
and Liu's symbolic factorization
algorithme nable s
ane #cie nt computation of
the minimale limination DAGs.
We also
de fine a
se t
dge s that must
be adde d to
the task-DAG in
orde r to
ge ne rate a minimal data-DAG that is valid as long as
partial pivoting with dynamic row and
columne xchange s is not
pe rforme d during
factorization. Finally,
we de scribe how
supple me nting this data-DAG
furthe r with a
small
se t
ofe xtrae dge s can
yie ld a
ne ar-minimal data-DAG that is
su#cie nt to
handle an arbitrary
numbe r of pivot
failure s and
the re sulting row and
columne xchange s
during
nume rical factorization. A pivot
failure occurs
the pivot
pre dicte d
by
the analysis
phase must
be alte re d during
nume rical factorization
be cause the nume rical
value of
the pivot is too small. By
me ans
ofe xpe rime nts on a
suite of
unsymme tric
sparse matrice s from
re al applications,
we show that computing
the final data-DAG
xtre me ly fast.
Furthe rmore , for
the matrice s in our
te st
suite , this
data-DAG has only a slightly
highe r
ofe dge s than
the task-DAG
constructe d
using
comple te transitive re duction.
The multifrontal
me thod [9, 14, 23] for
sparse matrix factorization usually
o#e rs
a significant
pe rformance advantage ove r
conve ntional factorization
sche me s by
pe r-
mittinge #cie nt utilization of
paralle lism and
me mory
hie rarchy. Du# and
Re id [14]
de scribe d a
symme tric-patte rn multifrontal algorithm for
unsymme tric
matrice s that
ge ne rate s
ane limination
tre e base d on
the symme tric
structure of
the union of
the structure s of A and
the transpose of A to
guide the nume rical factorization. This
algorithm works on
square frontal
matrice s
(se e se ction 4.1) and can incur a substantia

ove rhe ad for
ve ry
unsymme tric
matrice s
due to
unne ce ssary data
de pe nde ncie s
in
the e limination
tre e and
due toe xtra
ze ros in
the artificially
symme trize d frontal
matrice s. Davis and Du# [9] and
Hadfie ld [20]
introduce d an
unsymme tric-patte rn
multifrontal algorithm that
ove rcome s
the de ficie ncie s of a
symme tric-patte rn algo-
rithm. Our
powe rful symbolic
phase e nable s us to
use a much
more simplifie d and
e #cie nt
ve rsion of
the unsymme tric-patte rn multifrontal algorithm with partial pivoting

We de scribe the unsymme tric-patte rn multifrontal algorithm that is
use d in
WSMP
ande xpe rime ntally
compare it with
othe r
state -of-the -art
sparse unsymme tric
factorization
code s.
UNSYMMETRIC SPARSE MATRIX FACTORIZATION 531

Table
Test
ith their order (N), number of nonzeros (NNZ), and the application area of origin.
Number Matrix N NNZ Application
Finite element analysis
3 bayer01 57735 277774 Chemistry
4 bbmat 38744 1771722 Fluid dynamics
programming
6 e40r0000 17281 553956 Fluid dynamics
7 e40r5000 17281 553956 Fluid dynamics
simulation
9 epb3 84617 463625 Thermodynamics
engineering
14 mixtank 29957 1995041 Fluid dynamics
engineering
simulation
simulation
pre2 659033 5959282 Circuit simulation
19 raefsky3 21200 1488768 Fluid dynamics
22 tib 18510 145149 Circuit simulation
twotone 120750 1224224 Circuit simulation
wang3old 26064 177168 Circuit simulation
simulation
In

Table

1.1,
we introduce the suite of randomly
chose n
te st
matrice s that
we will
use ine xpe rime nts throughout this
pape r.
The table shows
the orde r
ofe ach matrix,
the numbe r of
nonze ros in it, and
the application
are a of
the origin of
the matrix. All
matrice s in our
te st
suite arise in
re al-life proble ms and
are in
the public domain.
The e xpe rime nts
re porte d in this
pape r
we re conducte d on an IBM RS6000 WH-2 with a
Gbyte s of RAM, 8
Mbyte s of
le ve l-2
cache , and 64
Kbyte s
of
le ve l-1
cache .
The organization of this
pape r is as follows.
ction 2
introduce s
the te rms,
conve ntions, and notations
use d in
the pape r. A symbolic factorization algorithm that
compute s
the structure of
the triangular factors and
minimale limination
structure s
is
de scribe d in
se ction 3. In
se ction 4,
we de scribe how to
compute ne ar-minimal
data-DAGs for
unsymme tric multifrontal factorization.
The nume rical factorization
algorithm is
discusse d in
de tail in
se ction 5.
We finish with concluding
re marks in
se ction 6.
The last
subse ction
ofe ach major
se ction
containse xpe rime ntal
re sults
pe rtaining to
the algorithms in that
se ction.
2. Terminology and conventions.
We assume that
the original n - n
sparse unsymme tric
coe #cie nt matrix is
irre ducible and cannot
be pe rmute d into a block-triangular
form. This is not a
se rious
re striction,
be cause a
ge ne ral matrix can first
be re duce d to a block-triangular form and
the n only
the irre ducible diagonal blocks
ne e d
to
be factore d [12].
We assume that
the coe #cie nt matrix A is
factore d into a
lowe r
triangular matrix L and an
uppe r triangular matrix U .
Multiple row and column
pe rmutations may
be applie d to A during various
stage s of
the solution
proce ss.
Howe ve r, for
the sake of clarity,
we will always
de note the coe #cie nt matrix by A and
the factors by L and U .
The state of
pe rmutation of A, L, and U will usually
be cle ar
from
the conte xt.
We de note the dire cte d graph
corre sponding to an n - n matrix M by
re graph may not always
be associate d
with
ane xplicitly
de fine d matrix.
Howe ve r,
the n
ane dge i#j # EM if
and only if m ij is a structural
nonze roe ntry in
the sparse matrix M .
The transpose of a matrix M is
re pre se nte d by M # . If i#j # EM ,
the n j#i # EM # , and
vice -ve rsa.
the se t of
indice s of
the columns in M that
have a structural
nonze roe ntry in row i. This is also
the se t of all
ve rtice s to which i has an outbound
e dge inG M . Similarly, Struct(M #,i ) is
the se t of
indice s of
the rows in M that
have a structural
nonze roe ntry in column i and is also
the se t of all
ve rtice s from which i
has an
inbounde dge inG M . A
dire cte d path from
node i to
node j in
the dire cte d
graphG M is
de note d by i#j.
The transitive re ductionG M O(VM , EM O ) of a graph
the graph with
the smalle st
dge s that has a
dire cte d path
i#j if and only
ifG M has a
dire cte d path i#j.
Since we are primarily
de aling with
the nonze ro
structure of
matrice s
rathe r than
the actual
value s,
we may also
loose ly
re fe r to M O as
the transitive re duction of M
O is a
transitive re duction
ofG M .
The le ading i - i submatrix of M is
de note d by M i and
the corre sponding graph and
its
transitive re duction
andG M O
re spe ctive ly.
The e dge s and paths in
some of
the graphs
use d in this
pape r
are labe le d. An
e dge in a
labe le d graph can
have one of
the thre e labe ls-L, U, or LU.
De pe nding on
its
labe l,
ane dge can
be an
L-e dge , a
U-e dge , or an
LU-e dge . L-, U-, and
LU-e dge s
from
ve rte x i to j
are de note d by i L
#j,
iU #j, and i
LU #j,
re spe ctive ly. An L-path from i
to j,
de note d by i L
#j, is a
dire cte d path containing only L- and
LU-e dge s. Similarly,
a U-path from i to j,
de note d by
iU #j, is a
dire cte d path containing only U- and
LU-e dge s. If an
L-e dge i L
#je xists in
the graph,
the
L-pare nt(i). Similarly, if
iU #je xists,
the
U-pare nt(i), and if i
LU #je xists,
the
LU-pare nt(i).
We de fine 1 a
supe rnode [q : r] as a maximal
se t of
conse cutive indice s {q, q
1, . , r} such that for all
and
matrice s L and U ,
we de fine
supe rnodal
matrice s L and U such
thate ach
supe rnode [q : r] in L
and U is
re pre se nte d by a
single row and column in L and U .
re m # n is
the total
numbe r of
supe rnode s.
Furthe rmore , if
and r < s,
the n g < h; that is,
the column and row
indice s in L and U maintain
the re lative orde r of
supe rnode s in L and U .
3. Computing a task-DAG and the structures of L and U .
Gilbe rt and
pre se nt an
unsymme tric symbolic factorization algorithm to
compute the structure s of
the factors L and U and
the ir
transitive re ductions L O and U O . Figuresummarize s
Gilbe rt and Liu's algorithm.
The algorithm
compute s
the structure
of L, U , and L O row by row and
compute s
the structure of U O by columns.
The total
time that
the algorithm shown in

Figure
spe nds in
ste p 1 is
bounde d
by flops(LU O ) [16], which is
the numbe r of
ope rations
re quire d to multiply
the sparse matrice s L and U O . Similarly,
the time spe nt in
ste p 3 is
bounde d by flops(UL O ).
The total computational cost of
ste ps 2 and 4 is
O| )). This is
be cause transitive re duction is
pe rforme d on n rows of U and columns of L, and
the ith
ste p
could
pote ntially
trave rse alle dge s
inG L O
andG U O
Ste ps 2 and 4 of
Gilbe rt and
Liu's algorithm
are much
more costly than
ste ps 1 and 3.
The cost of
the se ste ps
Other definitions of supernodes in the context of unsymmetric sparse factorization have been
used in the literature [11].
UNSYMMETRIC SPARSE MATRIX FACTORIZATION 533
to n do
1.
Compute
trave rsingG U O
and using
the fact that # j < i, j # Struct(L i,# ) if and only if # k # j such that k #
the re is a path k#j in U O
i-1 .
2.
Transitive ly
re duce Struct(L i,# )
usingG L O
ande xte nd it
toG L O
3.
Compute
4.
Transitive ly
re duce Struct(U #,i )
usingG U #O
ande xte nd it
toG U #O
end for
Fig . 3.1. Gilbert and Liu's unsymmetric symbolic factorization algorithm [16].
to n do
1.
Transitive ly
re duce Struct(L i,# )
usingG L O
ande xte nd it
toG L O
2.
Compute
3.
Transitive ly
re duce Struct(U #,i )
usingG U #O
ande xte nd it
toG U #O
4.
Compute
end for
Fig . 3.2. A modified symbolic factorization algorithm.
has
prompte d
re se arche rs to
se e k
alte rnative s, such as computing fast but
incomple te transitive re duction [9, 15].
The use of such
alte rnative s
toG L O
andG U O with
more e dge s
thanG L O
andG U O ,
re spe ctive ly, can
incre ase the cost of
ste ps 1 and 3, as
we ll
as that of
nume rical factorization.
3.1. A modification to Gilbert and Liu's algorithm.
We now
de scribe a
re lative ly
simple modification to
the algorithm shown in

Figure

3.1.
We start by
splitting
the original
coe #cie nt matrix into a
lowe r triangular part
store d by columns
and an
uppe r triangular part
store d by rows. In our
modifie d symbolic factorization
algorithm,
we compute the structure of L by
the columns
(i.e ., L # by rows) and that
of U by
the rows. This is
achie ve d by simply
re formulating
the algorithm shown in

Figure

3.1 to
pe rform only
ste ps 2 and 3, but
twice fore ach i on two
se ts of
ide ntical
data
structure s-one corre sponding to L # and
the othe r
corre sponding to U .
The modifie d algorithm is shown in

Figure

3.2.
Note that in
the algorithm of

Figure

3.2,
ste ps 3 and 4
are ide ntical to
ste ps
1 and 2,
re spe ctive ly.
The first two
ste ps
compute the ith rows of L O and U and
the last two
ste ps
compute the ith columns of U O and L. An actual
code of this
algorithm can
use the same pair of
routine s with
di#e re nt
argume nts to
imple me nt
all four
ste ps.
The re duction in
the size of
the code by half,
howe ve r, is a
se condary
be ne fit of
the modifie d algorithm.
The primary
advantage of this
sche me is that it
allows
imme diate de te ction of
supe rnode s during symbolic factorization. This, as
we shalle xplain in
se ction 3.2, allows us to avoid computing and
storingG L O
andG U #O
e xplicitly.
Inste ad,
we can work only with
the ir
supe rnodal
counte rpartsG L O and
G U #O .
3.2. Use of supernodes to speed up transitive reduction. Most
mode rn
sparse factorization
code s
re ly
he avily on
toe #cie ntly
utilize me mory
hie rarchie s and
paralle lism in
the hardware .
are so crucial to high
pe r-
formance in
sparse matrix factorization that
the crite rion for
the inclusion of rows
and columns in
the same supe rnode is
ofte n
re laxe d [7] to
incre ase the size of
the supe rnode s.
Conse cutive rows and columns with
ne arly
the same but not
ide ntical
structure s
are ofte n
include d in
the same supe rnode , and artificial
nonze roe ntrie s with
a
nume rical
value of 0
are adde d to maintain
ide ntical row and column
structure s for
all
me mbe rs of a
supe rnode .
The rationale is that
the slight
incre ase in
the numbe r
of
nonze ros and floating-point
ope rations
involve d in
the factorization is
more than
compe nsate d for by a
highe r factorization
spe e d.
WSMP's LU factorization algorithm also works on
the re laxe d
ge ne rate
d by its symbolic factorization. In
the symbolic factorization algorithm, as soon as
are compute d in
the ith
ite ration of
the oute r loop,
the y
can
be compare d with Struct(L #,i-1 ) and Struct(U i-1,# ) to
de te rmine if
the y
be long
to
the curre nt
supe rnode . A
ne w row-column pair is
adde d to
the curre nt
supe rnode if its
structure
ide ntical or
ne arly
ide ntical to
the pre vious row-column pair.
If
the ith row-column pair fails to
me e t
the crite rion for
me mbe rship into
the curre nt
supe rnode ,
the n a
ne w
supe rnode is
starte d at i.
The use of
supe rnode s allows us to significantly
re duce the cost of computing
the transitive re ductions. In
ste p 1 of
the algorithm shown in

Figure

3.2,
inste ad of
transitive ly
re ducing
the e ntire Struct(L i,# ),
we re duce only
the se t {h
re
Ste p 3 is
tre ate d similarly. As a
re sult of working only
with
supe rnode s,
the uppe r bound on
the cost of computing
the transitive re duction
de cre ase s from
O| )) to
O| )). This is
be cause only
the supe rnodal
DAGsG L O
andG U O
are se arche d
duringe ach of
the n
transitive re duction
ste ps. Strict
supe rnodal
andG U O would
have at
le ast
fe we re dge s
thanG L O
andG U O ,
re m is
the numbe r of
supe rnode s.
The re ason
is that U O and L #O do not contain
anye dge s i#j,
re
supe rnode .
The use of
re laxe d
re duce s
the numbe r
ofe dge s
e ve n
furthe r
be cause some pote ntiale dge s of
the form i#j,
re
be e liminate d from
the task-DAG
node s i and
are artificially
me rge d.
3.3. Task-DAGs for LU factorization. In this
pape r,
we will
re fe r to two
type s of task-DAGs: a
conve ntional DAG
de note d by
and a
supe rnodal DAG
de note d by T S . Each
ve rte x of
the conve ntional task-DAG
re fe rs to
the task of
computing a
single row of U and
the corre sponding column of L. On
the othe r hand,
a
ve rte x of
the supe rnodal task-DAG
corre sponds to a
se t of row-column pairs that
constitute a
supe rnode . Although, in a practical
imple me ntation,
we always work
with
supe rnodal DAGs,
we will
ofte n
use conve ntional task- and data-DAGs in
the re mainde r of
the pape r to
the e xposition
simple . All
re sults and
de scriptions
pre se nte d in
te rms of
the conve ntional DAGs map naturally to
the supe rnodal
case .
We first show how to
compute TC in
te rms of
the conve ntional
structure s L #O
and U O .
The transpose matrix L # is
use d to
indicate that for all i#j
Theorem 3.1.
is a task-DAG for LU factorization if its
ve rte x
se
{1, 2, . , n} and
itse dge -se
Proof. To
prove that
is a task-DAG,
we show that E T C is
su#cie nt to
re p-
re se nt a
prope r
orde ring of
the ne limination tasks
de note d by V T C . Struct(L #,i ) can
contribute to Struct(L #,j ) only if i # Struct(U #,j ), and if this is
the case ,
the n
the symbolic factorization algorithm of

Figure

3.2e nsure s that U O
containse ithe r i#j
UNSYMMETRIC SPARSE MATRIX FACTORIZATION 535

Table
Comparison of conventional symbolic factorization (due to Gilbert and Liu
[16])w ith supernodal
| is the size of the largest diagonal block in the matrix
onw hich symbolic
factorization is performed; Nsup is the number of supernodes; t C and t S are the times in seconds
of the
tw are the number of edges in the
task-DAGs produced by the
tw algorithms.
Matrix |V | Nsup
Conventional
Supernodal
|t Ct S
bbmat 38744 4877 10. 41260 1.7 6077 7.9 6.8 5.9
e40r5000 17281 2755 .60 19891 .16 3182 6.3 6.3 3.8
fidap011 16614 1262 2.3 16613 .42 1261 13. 13. 5.5
mil053 530238 166155 15. 530237 4.5 166154 3.2 3.2 3.3
mixtank 29957 2984 7.8 30949 1.2 3203 10. 9.7 6.5
nasasrb 54870 3808 4.9 54869 .97 3807 14. 14. 5.1
pre2 629628 243693 30. 765210 6.4 317216 2.6 2.4 4.7
raefsky3 21200 1282 2.1 21199 .41 1281 17. 17. 5.1
raefsky4 19779 1359 2.9 19778 .50 1358 15. 15. 5.8
tib 17583 7823 .11 22904 .07 10060 2.2 2.3 1.6
twotone 105740 34304 2.6 126656 .91 44856 3.1 2.8 2.9
wang3old 26064 8451 3.1 26063 .54 8450 3.1 3.1 5.7
wang4 26068 8254 3.0 26067 .53 8253 3.2 3.2 5.7
or i#j.
The same is
true for Struct(U i,# ), Struct(U j,# ), and L #O .
The re fore ,e ve ry
row-column pair i that
update s row and column j must
be e liminate d
be fore j.
The ore m 3.1 can
be e asilye xte nde d to
the supe rnodal
case .
The supe rnodal
task-DAG T S is
de fine d by a
ve rte x
se
ane dge se
re m is
the numbe r of
supe rnode s.
3.4. Experimental results. In

Table

3.1,
we compare Gilbe rt and Liu's symbolic
factorization algorithm [16] with
the supe rnodal symbolic factorization algorithm
de scribe d in
se ction 3.2.
We re port
the ir CPU
time s
tC and t S ,
re spe ctive ly, and
the numbe re dge s in task DAGs
and T S
ge ne rate d by
the m.
The last column of

Table

3.1 shows
the factor by which
the supe rnodal symbolic
factorization is
faste r than
the conve ntional algorithm.
The table also shows
ave rage supe rnode size (n/m) and
the ratio
dge s in
and T S
fore ach matrix.
The se two
ratios
are close ly
re late d.
The ratio of
tC and t S
be ars
some corre lation to
the ratio
dge s in
and T S , but
the actual ratio is matrix
de pe nde nt.
Note that only
the time of
transitive re duction
ste ps 1 and 3 of
the algorithm in

Figure

3.2 is
re duce d by
the use of
supe rnode s;
the time of computing
the structure s of L and U in
ste ps 2 andre mains mostly
unchange d
(othe r than
some re duction in
the numbe r of
structure s
me rge d
due to
supe rnode re laxation).
The re fore ,
the actual
re duction
achie ve d in
the symbolic factorization
time de pe nds on
the re lative amounts of
time spe nt in
536 ANSHUL
transitive re duction and computing L and U
structure s.
More ove r,

Table
re ports
the numbe r
dge s in
the task-DAGs, not
the numbe r
dge s in
the actual
lowe r
and
uppe r triangular
transitive ly
re duce d graphs that
are trave rse d during symbolic
factorization.
Re call that
the e dge -se t of a task-DAG is
the union of
the e dge -se ts
of
the corre sponding
lowe r and
uppe r triangular
transitive ly
re duce d graphs.
The amount of structural
symme try in
the matrix
a#e cts
the numbe r of
commone dge s
be twe e n
the uppe r and
lowe r
transitive ly
re duce d graphs, which in turn
de te rmine s
the actual
dge s in
the task-DAG.
Eise nstat and Liu [15]
pre se nt an
alte rnative to
comple te transitive re duction
to
re duce the cost of this
ste p in
sparse unsymme tric symbolic factorization.
The y
propose e xploiting structural
symme try in
the matrix to
compute partial
transitive re ductions. Although
the y
pre se nte xpe rime ntal
re sults on a
di#e re nt
se t of much
smalle r
matrice s, it
appe ars that
the use of
supe rnode s as
propose d in
se ction 3.2
can
achie ve much
highe r
spe e dups in symbolic factorization
while computinge xact
transitive re ductions than
the partial
transitive re duction
sche me propose d in [15].
Howe ve r,
Eise nstat and Liu's algorithm too can
be spe d up by
the use of
supe rnode s.
A
supe rnodal
ve rsion of this algorithm has
be e n
imple me nte d in
the Supe rLU dist [21]
sparse solve r
package .
We compare d our symbolic factorization
time with that of
Supe rLU dist and found
the latte r to
be slowe r by about 25%
ove rall on our
te st
suite . This could
be partly
due to
imple me ntation
di#e re nce s and partly
due to
the fact that
while Eise nstat and Liu's algorithm
save s
time in
the transitive re duction
computation, it
spe ndse xtra
time in
me rging
structure s
due to
re dundante dge s in
the DAG. It
appe ars that
the use of
supe rnode s in
Gilbe rt and Liu's algorithm can
spe e d
up its
transitive re ductione nough for it to match or
outpe rforme ve n a
supe rnodal
ve rsion of
Eise nstat and Liu's algorithm
ine xe cution
time .
4. Data-DAGs for unsymmetric multifrontal LU factorization.
The original
multifrontal algorithm [14, 23] was
de scribe d in
the conte xt of a
symme tric-
patte rn
coe #cie nt matrix but has
be e n
applie d to
matrice s with
unsymme tric patte
rns by introducing
ze ro-value de ntrie s at
appropriate locations to
conve rt
the original
matrix into
one with
the patte rn of A+A # [14, 2, 4]. This can
cause a substantial
ove rhe ad for
ve ry
unsymme tric
matrice s
due to
the e xtra computation
pe rforme d on
the introduce de ntrie s and
the re sulting fill-in. Davis and Du# [9] and
Hadfie ld [20]
introduce d an
unsymme tric-patte rn multifrontal algorithm to
ove rcome this short-
coming. In this
se ction,
we de ve lop
ne ar-minimal data-DAGs for
the unsymme tric
multifrontal algorithm-an
aspe ct of
unsymme tric multifrontal factorization that has
not
be e n
we ll
inve stigate d in
pre vious works. As
we shall show in
se ction 5,
the availability of a
ne ar-minimal data-DAG aids in
the e #cie nt
imple me ntation of
the nume rical factorization
phase . It would also
he lp
minimize the synchronization and
communication
ove rhe ads in a
paralle l
imple me ntation.
4.1. Outline of the symmetric multifrontal algorithm.
The symme tric-
patte rn multifrontal algorithm is
guide d by an
asse mbly
ore limination
tre e [22, 23,
19], which
se rve s as both
the task- and
data-de pe nde ncy graphs for
the factorization
proce ss.
The data
associate d
withe ach
supe rnode of
the e limination
tre e is a
square frontal matrix. A frontal matrix F g
associate d with a
supe rnode
de nse matrix
whose dime nsions
are e qual
)| or| Struct(U q,#
)| .
The contiguous
local row and column
indice s in
the de nse frontal matrix
corre spond to noncontiguous
global
indice s of
the matrix L
Eache ntry in a frontal matrix
corre sponds to
a structural
nonze roe ntry in
the global matrix.
Afte r a frontal matrix F g is fully
asse mble d or
populate d,
the le ading r - q columns
corre sponding to
UNSYMMETRIC SPARSE MATRIX FACTORIZATION 537
the supe rnode (also known as
the pivot block)
are factore d and
be come parts of
the factors U and L,
re spe ctive ly.
The re maining trailing part of
the frontal matrix is now
calle d
the update or
the contribution matrix,
de note d by C g .
The contribution matrix
corre sponding to a
supe rnode is
asse mble d
comple te ly into
the frontal matrix of its
only
pare nt
supe rnode and is
ne ve r
acce sse d again. This is
be cause if
the pare nt of
supe rnode
the e limination
tre e ,
the n Struct(L #,r
The same is
true for columns of U
due to
symme try.
In a
re cursive formulation of
the symme tric-patte rn multifrontal algorithm,
the task
corre sponding to a
supe rnode first
comple te s
ide ntical subtasks
fore ach of its
childre n in
the e limination
tre e ,
the n
asse mble s
the ir contribution
matrice s into its
frontal matrix, and finally
pe rforms
the partial factorization on
the frontal matrix.
Calling a
re cursive proce dure to
pe rform
the task
de scribe d
above on
the root su-
pe rnode of
the e limination
tre e comple te s
the factorization of a
sparse matrix with a
symme tric
structure .
4.2. Outline of the unsymmetric multifrontal algorithm.
The ove rall
structure of an
unsymme tric-patte rn multifrontal algorithm is similar to its
symme t-
ric
counte rpart and can
be e xpre sse d in
the form of a
re cursive proce dure starting at
the root
(the supe rnode with no
outgoinge dge s) of
the task-DAG.
Howe ve r,
the re are two major
di#e re nce s.
The first
di#e re nce is in
the control-flow. In
the unsymme tric
multifrontal algorithm,
be fore starting a subtask for a child,
the task
corre sponding
to
the pare nt
supe rnode must
che ck to
se e if
the child
supe rnode has
alre ady
be e n
proce sse d by
anothe r
pare nt. Only
the first
pare nt to
re ach a child actually
pe rforms
the re cursive computation starting at that child.
The se cond
di#e re nce is in
the data-
flow, or
the way contribution
matrice s
are asse mble d into frontal
matrice s. This is
e xplaine d
be low in
gre ate r
de tail.
Re call that
the e dge -se t E T C of
the task-DAG
is
the union of
the e dge -se ts
O of
the transitive re ductions of L # and U ,
re spe ctive ly.
We now assign
labe ls to
the e dge s in
.
The e dge s
contribute d to E T C
sole ly by E L #O
are labe le d as
L-e dge s.
dge s
contribute d to E T C
sole ly by E U O
are labe le d as
U-e dge s.
The third
type of
labe l,
the LU-labe l, is
assigne d to
the e dge s that
be long to
the inte rse ction E L #O and E U O . Finally, an
L-e dge i L
#j is
conve rte d to an
LU-e dge i
LU #j
if
the re is a U-path
iU #j in
, and a
U-e dge iU #j is
conve rte d to i
LU #j if
the re is an
L-path i L
#j in
.
The e dge s of
the supe rnodal task-DAG T S
are de fine d similarly.
Unlike the symme tric multifrontal algorithm,
the frontal and contribution matrice
s in
the unsymme tric multifrontal algorithm
are , in
ge ne ral,
re ctangular
rathe r
than
square .
Furthe rmore , a contribution matrix in
the unsymme tric multifrontal
algorithm can
pote ntially
be asse mble d into
more than
one frontal matrix
be cause a
supe rnode in
the data-DAG can
have more than
one pare nt. As
de scribe d in [20],
the asse mbly of contribution
matrice s into
the pare nt frontal
matrice s in
the unsymme tric
multifrontal algorithm
proce e ds as follows.
#h
be an
L-e dge in
the data-DAG,
re
have an
inde x i in common,
the n
alle le me nts of row i
of U in C g can
pote ntially
be asse mble d into F h . Similarly, if
gU #h is a
U-e dge and
have an
inde x i in common,
the n
alle le me nts of column
i of L in C g can
pote ntially
be asse mble d into F h . Finally, if g
LU #h is an
LU-e dge ,
the n
the e ntire trailing submatrix of C g with global row and column
indice s
gre ate r
than
ore qual to s can
be asse mble d into F h .
Ce rtaine ntrie s of C g may
have pote ntial
de stinations in
the frontal
matrice s of
Multifrontal factorization guided by the task-DAG
Matrix
LU
U
LU
LU
LU U
intended destination
task-DAG
U
Fig . 4.1. An example of the inability of a task-DAG to guide complete assembly of all contribution
matrices in the unsymmetric multifrontal algorithm. An 'X' denotes a nonzero in the coe#cient
matrix and a '+' denotes a nonzero created due to fill-in.
more than
one pare nt of
ge ve n if
the data-DAG contains no
unne ce ssarye dge s. This
is
be cause C g can
have common rows (columns) with
the frontal
matrice s of
more than
one among g's LU- and
L-pare nts
(U-pare nts).
The unsymme tric multifrontal
algorithm
muste nsure that
anye ntry of a contribution matrix is not
use d to
update more than
one frontal matrix. Additionally, a
corre ct data-DAG must
have su#cie nt
outgoinge dge s from all
supe rnode s so
thate ache ntry of a contribution matrix has a
pote ntial
de stination in at
le ast
one frontal matrix.
4.3. Inadequacy of task-DAG for unsymmetric multifrontal algorithm.
By
me ans of a
smalle xample in

Figure

4.1,
we show that if
the task-DAG
de fine d
in
se ction 3.3 is
use d as a data-DAG,
the n all contribution
matrice s may not
be fully
absorbe d into
the ir
pare nt frontal
matrice s.
The figure shows a
sparse matrix
with factorization fill-in,
the transitive ly
re duce d DAGs L #O and U O , and
the task-
DAG with
itse dge s
labe le d as
de scribe d in
se ction 4.2. For
the sake of
clarity,e ach
supe rnode is
chose n to
be of
size 1.
The figure shows all frontal and contribution
(shade d portions)
matrice s and
the flow of data from
the contribution to frontal
matrice s along
the e dge s of
the task-DAG.
Note that
alle dge s may not
le ad to a data
LU #5. It
asily
se e n that
the U-e dge 1U #4, which is
abse nt from
the task-DAG
(be cause it is
re move d
while transitive ly
re ducing U to U O ), is
ne ce ssary
for
the comple te asse mbly of C 1 .
4.4. A data-DAG for a predefined pivot sequence. Having shown that
the minimal task-DAG cannot
se rve as a data-DAG for
unsymme tric multifrontal
we now
de fine a data-DAG that is
su#cie nt for
the prope r
asse mbly
of all contribution
matrice s, as long as rows and columns
are note xchange d among
di#e re nt
supe rnode s for pivoting.
We will
use D N to
de note such a DAG,
whe re the supe rscript N stands for "no pivoting." A data-DAG D P that can
accommodate
UNSYMMETRIC SPARSE MATRIX FACTORIZATION 539
pivoting will
be de scribe d in
se ction 4.5.
Theorem 4.1. If a column
s all of
the following
conditions,
the n a
U-e dge i U
#j is
ne ce ssary for C i to
be comple te ly
asse mble d into its
pare nts' frontal
matrice s:
1.
The LU-pare nt of i, if
ite xists, is
gre ate r than j.
2.
None of i's
U-pare nts
are in Struct(U #,j ).
3.
The re e xists a k # Struct(L #,i ) such that k > j.
The transpose of this
the ore m can
be state d similarly.
Proof.
The contribution matrix C i has a column that
contribute s to L #,j ,
be cause ,
at
the le ast,
the re is
ane le me nt
corre sponding to L k,j in C i . At
the same time ,
none of i's
U-pare nts' frontal
matrice s
have column j, so
the y cannot absorb L #,j from C i .
Since the LU-pare nt of i is
gre ate r than j, it too cannot absorb L #,j from C i .
The addition of
iU #j
make s it
possible for C i to
contribute L #,j to F j .
The transpose case can
be prove n similarly.
The ore m 4.1
capture s
the situation
illustrate d in

Figure

4.1 for 4, and
pre scribe s
the addition of
toe nsure comple te asse mbly of C 1 .
Theorem 4.2. If D N is a DAG
forme d by adding all
possible e dge s to
according to
The ore m 4.1,
provide d that
the se e dge s don't
alre adye xist,
the n D N is a
data-de pe nde ncy DAG for
the unsymme tric multifrontal algorithm without pivoting.
Proof. To show that D N is a data-DAG,
we must show that
itse dge -se t is
su#cie nt for
the comple te absorption of all contribution
matrice s into
the ir
pare nt
frontal
matrice s.
We prove this by contradiction.
Without loss of
ge ne rality,
assume that
ane le me nt
corre sponding to L k,j in C i is
not
asse mble d.
Note that i < j < k. If L k,j is in C i ,
the
ithe r
the re is a U-path
iU #j in
.
If
the n
alle ntrie s with row
indice s
gre ate r than
ore qual to j in column
of C i will
be absorbe d by F j , and
the se e ntrie s
include the one corre sponding to
L k,j . If
the n a U-path
iU #je xists in
the re are two
possibilitie s:
e ithe r
LU-pare
LU-pare nt(i) > j.
LU-pare nt(i). If l # j,
the n
the e ntire trailing submatrix of C i with row and column
indice s
gre ate r than l, including
be asse mble d into F l . If l > j,
the n
conside r two
furthe r
possibilitie s:e ithe r
one of i's
U-pare nts is in Struct(U #,j ) or is not. If
one is,
the n its frontal matrix
will absorb column j from C i . If
none of i's
U-pare nts is in Struct(U #,j ),
the n all
conditions for
the applicability of
The ore m 4.1
are satisfie d.
The re fore ,
iU #j would
have be e n
adde d to D N and would
have cause d
the e ntry
corre sponding to L k,j in
C i to
be absorbe d into F j . Thus, it is not
possible for
the e ntry
corre sponding to
L k,j to
be le ft
unasse mble d in any C i . Similarly, it can
be shown that
the e ntry
corre sponding to any U j,k cannot
be le ft
unasse mble d in any C i .
Having shown that
the e dge -se t of D N is
su#cie nt for
unsymme tric multifrontal
factorization without pivoting,
we now show that not
alle dge s that D N
inhe rits from
may
be ne ce ssary if pivoting is not
pe rforme d during factorization.
Theorem 4.3. For LU factorization without pivoting,
ane dge i U
#j (i L
in
is
re dundant if
the maximum
inde x in Struct(L #,i ) (Struct(U i,# )) is
smalle r
than j.
Proof.
Re call that Struct(L #,j
the maximum
inde x in Struct(L #,i ) is
smalle r than j,
the n
doe s not
contribute to Struct(L #,j ).
The proof for
L #O and Struct(U i,# ) is similar.
Note that
The ore m 4.3 is valid only if row and
columne xchange s
are not
pe r-
forme d during LU factorization.
Othe rwise , additional fill-in
cause d by pivoting could
cre ate an
inde x
gre ate r than
ore qual to j in Struct(L #,i ) or Struct(U i,#
),e ve n if it is
not
pre dicte d by
the symbolic factorization on
the original
pe rmutation of
the matrix.
The re fore ,
alle dge s in
could
pote ntially
be use d.
Supe rnodal
ve rsions of
The ore ms 4.1-4.3 for T S can
be prove n similarly. To
summarize the re sults of this
subse ction,
we have shown how to construct a data-
DAG for
unsymme tric multifrontal factorization without pivoting from a task-DAG
and
we have shown that although
the task-DAG is
de rive d from
the strict
transitive re ductions of L # and U (or L # and U), it may still pass
one dge s to
the data-DAG
that
are re dundant if pivoting is not
pe rforme d during factorization.
The re fore ,
the data-DAG is not minimal.
Howe ve r, if pivoting is
pe rforme d,
the n
pote ntially all
the e dge s could
ge t
use d.
4.5. Supplementing the data-DAG for dynamic pivoting.
We will now
show that
the e dge -se t of data-DAG D N
constructe d in
se ction 4.4 may not
be su#-
cie nt if pivoting is
pe rforme d during factorization.
We also discuss how to
supple me nt
to
ge ne rate a data-DAG D P
whose e dge -se t is
su#cie nt to
handle any amount of
pivoting.
We start with an
ove rvie w of
the pivoting
me thodology in
the unsymme tric
multifrontal algorithm, which has
be e n
de scribe d in
de tail in [20].
If a
diagonale le me nt A i,i (q # i # r) in a
supe rnode [q : r] fails to
me e t
the pivoting
crite rion,
the n first an
atte mpt is
made toe xchange row and column i with
a row j and a column k such that i <
the pivoting
crite rion. Such
intrasupe rnode pivoting has
noe #e ct on
the structure of
the factors
and factorization can
continue as usual.
Howe ve r, it may not always
be possible to
find a
suitable row-column pair within a
supe rnode 's pivot block to satisfy
the pivoting
crite rion. In this situation,
inte rsupe rnode pivoting is
ne ce ssary. If
the LU-pare nt of
the data-DAG and a
suitable ith pivot cannot
be found
within
the pivot block of F g ,
the n all row-column pairs from i to r
are symme trically
pe rmute d to
ne w locations from s - (r
Thus,e #e ctive ly,
shrinks to [q
the supe rnode [s :
t]e xpands to [s - (r
As a
side e #e ct of this pivoting,
the re is additional fill-in in all
the ance stors of g in
the data-DAG that
are smalle r than h. In particular,
the columns of L of all of g's
U-ance stors
smalle r than h
ge te xtra row
indice
the rows of U of all of g's
L-ance stors
smalle r than h
ge te xtra column
indice
failure in
supe rnode h is
handle d similarly in a
re cursive manne r.
In D N ,
whose construction is
de scribe d in
se ction 4.4, all
supe rnode s may not
have an
LU-pare nt to support
the symme tric pivoting
me thod
de scribe d
above .
The re fore ,
as
the first
ste p towards
de riving D P from D N ,
we alte r
the e dge -se t of
the latte r
as follows.
Fore ach g from 1 to m
re m is
the total
numbe r of
supe rnode s),
the smalle st
supe rnode h to which both g L
#h and
gU #he xist is
de signate d as
the LU-pare nt of g; that is, if
ane dge g#h
doe s
note xist,
the n an
LU-e dge g
LU #h is
adde d to
the data-DAG, or if an L- or a
U-e dge g#he xists,
the n it is
conve rte d to an
LU-e dge .
The n,
alle dge s g#k such that k > h
are de le te d. If
the original matrix is
not
re ducible to a block-triangular form,
the n
afte r this
modification,e ach
supe rnode othe r than
the root
supe rnode has an
LU-pare nt to
accommodate row-column pairs
that fail to satisfy
the pivoting
crite rion in
the ir original locations [20]. It
asily
se e n that this modification has
noe #e ct on
The ore ms 4.1-4.3
be cause g L
#h
(gU #h) is
in
the modifie d D N only if g L
#h
(gU #h) is in
the original D N as
de fine d in
se ction 4.4.
UNSYMMETRIC SPARSE MATRIX FACTORIZATION 541
old
indices
Original matrix
Multifrontal factorization without any pivot failure
FLU
Factorization with failure of pivot 1
Handling failure of pivot 1
new
unassembled intended destination
data-DAG13LU
U
LU
LU
Fig . 4.2. An example factorization to
show how the failure of pivot 1 is handled by a symmetric
permutation of
row and column 1 to merge
themw ith their LU-parent supernode, 4. An 'X' denotes
a nonzero in the coe#cient matrix and a '+' denotes a fill-in. The circled `X' and '+' are created
due to pivoting. A '0' denotes a fill-in predicted by the original symbolic factorization that has a
value of zero due to pivoting-related movement of
row s and columns. The figure also
show s that the
absence of 2 L
#4 leaves the entry U 1,5 unassembled from C 2 .

Figure

4.2 shows how
the failure of pivot row and column 1 is
atte mpte d in
the unsymme tric multifrontal factorization of a small 5 -
5e xample matrix. Row-column
1 is
symme trically
pe rmute d to a
ne w location
adjace nt to 1's
LU-pare nt 4 in
the data-
DAG. This
re sults in an addition of row
inde x 1 to 1's
U-pare nt 2 and an addition of
column
inde x 1 to 1's
L-pare nt 3. Additionally,
afte r moving to
the ir
ne w locations,
row 1 in U and column 1 in L
ge t fill-in in column and row positions
re row 4 in
U and column 4 in L
have nonze ros
(i.e ., U 1,5 , L 4,1 , and L 5,1 ).

Figure

4.2 also shows
that
afte r pivoting,
the ne w row 1 of C 2 cannot
be fully
asse mble d in
the abse nce of an
L-e dge 2 L
#4.
arly, in addition to adding
dge s as
de scribe de arlie r, D N
re quire s
furthe r modifications in
orde r to
se rve as a data-DAG for
unsymme tric multifrontal
algorithm with dynamic pivoting.

Figure
give s
anothe re xample of a factorization
whe re D N is
unable to
guide a
comple te asse mbly in
the e ve nt of a pivot
failure .
Note that
thise xample satisfie s
the first two conditions of
The ore m 4.1.
Howe ve r,
since it
doe s not satisfy condition 3, no
dge s
are adde d and
the abse nce of
pre clude s a
comple te asse mbly of C 2 into its
pare nts' frontal
matrice s
We now
state and
prove a
the ore m that
pre scribe s a modification of D N to
pre ve nt
the situation
illustrate d in

Figure

4.3.
Theorem 4.4. If a column
s all of
the following
conditions,
the n a
U-e dge i U
#j is
ne ce ssary for C i to
be comple te ly
asse mble d into its
pare nts' frontal
matrice s in
the e ve nt of
failure of pivot k.
1.
The LU-pare nt of i is
gre ate r than j.
2.
None of i's
U-pare nts
are in Struct(U #,j ).
3. A
xists such that
the re is a U-path k U
#i in
LU-pare nt(k) > j.
The transpose of this
the ore m can
be state d similarly.
Proof.
Note that
The ore m 4.4 is
ve ry similar to
The ore m 4.1.
The only
di#e re nce is condition 3. If pivot k fails,
the n it will add a row in Struct(L #,i ) that
corre sponds
to
LU-pare nt(k) -1, which is
the ne w location of k and is
gre ate r than j - 1,
the ne w
inde x for j. Thus,
the failure of pivot k transforms condition 3 of
The ore m
4.4 into condition 3 for
the applicability of
The ore m 4.1, which has
alre ady
be e n
prove d.
The ore m 4.4
state s
thate ve n if Struct(L #,i )
doe s not
have any
inde x
gre ate r than
j but all
othe r conditions for
the applicability of
The ore m 4.1
are satisfie d and
iU #j
is not
pre se nt in
the DAG,
the n pivoting may
re sult in
incomple te asse mbly
unle ss
thise dge is
adde d. This is
be cause pivoting can
cre ate a
nonze roe ntry L k,i such that
j. This is what
happe ns in
the e xample shown in

Figure

4.3 for
the original
indice s. Pivoting
change s i, j, and k to 1, 7, and 8,
re spe ctive ly.
In light of
The ore m 4.4,
we introduce anothe r modification to D N .
Inste ad of using
The ore m 4.1 strictly to
de rive D N from
we omit
che cking for condition 3 and
de rive D N by adding all
those e dge s to
practice ) that satisfy conditions 1
and 2.
Now, by
me ans of
The ore m 4.5,
we will show that
the data-DAG D N
,e ve n
afte r
the modifications
de scribe d
above , is not
su#cie nt
toe nsure comple te asse mbly of
all contribution
matrice s in
the e ve nt of
inte rsupe rnode pivoting.
The re ade r can
ve rify that

Figure
the transpose case of
The ore m 4.5 for
5. Finally,
The ore m 4.6 will show that
supple me nting
the data-
DAG with
additionale dge s
pre scribe d by
The ore m 4.5
make s it
su#cie nt to
handle all contribution
matrice s in
the face of
inte rsupe rnode pivoting. As
we dide arlie r in
this
pape r, for
the sake of clarity and simplicity,
we will
state and
prove The ore ms 4.5
and 4.6 in
the conte xt of
conve ntional DAGs with
single -node supe rnode s.
The re sults
naturallye xte nd to
supe rnodal DAGs.
Theorem 4.5. If h is
the LU-pare nt of j and all of
the following conditions hold,
the n a
U-e dge i U
#h is
ne ce ssary for C i to
be comple te ly
asse mble d into its
pare nts'
frontal
matrice s in
the e ve nt that j fails to
me e t
the pivot
crite rion in its original
location.
1.
The re e xists an L-path j L
#i such that i < h and
LU-pare nt(i) > h.
2.
None of i's
U-pare nts
are in Struct(L #,j ).
3.
Eithe r # k # Struct(L #,i ) such that k > h, or
the re is a U-path k U
LU-pare
The transpose case can
be state d similarly.
UNSYMMETRIC SPARSE MATRIX FACTORIZATION 543
indices
indices
new
F 9
F 486423748
9 LU
UU
U
U
U
U
U
Original matrix data-DAG9 X
unassembled5
Handling failure of pivot 1
intended destination
LU
LU
LU
LU
LU
LU
Fig . 4.3. An example factorization to
show that the edges in D N are not su#cient to assemble
its parents' frontal matrices in the event of the failure of pivot 1. The convention for
representing di#erent types of structural nonzeros is the same as in Figure 4.2.
Proof. If pivot j fails,
the n, along with
the othe r
faile d
LU-childre n of h, it
occupie s a
ne w position just
be fore h.
Since the re is an L-path j L
#i, column j is
adde d to C i
afte r
the failure of pivot j; that is, in
the ne w matrix
afte r pivoting, j #
We know that
the LU-pare nt of i is
gre ate r than
the ne w j,
be cause LU-pare nt(i) > h.
Since none of i's
U-pare nts
we re in
the old Struct(L #,j ),
the y
are not in
the ne w Struct(U #,j
)e ithe r. Thus
the first two conditions for
the applicability
of
The ore ms 4.1 and 4.4
are satisfie d. Condition 3 of
The ore m 4.5
quivale nt to
condition 3 of
The ore ms 4.1 and 4.4.
The re fore , a
U-e dge iU #j is
ne e de d for
prope r
multifrontal factorization of
the ne w matrix
afte r
pe rmuting j to its
ne w location.
Since , in its
ne w location, j is
me rge d with h into a common
supe rnode , a
U-e dge iU #h in
the original matrix would
have su#ce d.
The proof of
the transpose case is
similar.
Theorem 4.6. If D P is a DAG
forme d by adding all
possible e dge s according to
The ore m 4.5 to D N ,
the n D P is an
ade quate data-DAG for
unsymme tric multifrontal
factorization with
pote ntially
unlimite d
inte rsupe rnode pivoting.
Proof.
We prove this by showing that with D P , it is not
possible for
anye le me nt of
a contribution matrix C i to
re main
unasse mble d. Without loss of
ge ne rality,
conside r
ane le me nt
corre sponding to L k,j in C i . If L k,j is in C i ,
the ne ithe r k # Struct(L #,i )
the original L and U
pre dicte d by symbolic factorization, or
row k or column j or both
we re adde d to C i
due to pivoting. If row k and column j
are parts of
the original
structure of C i ,
the n
The ore m 4.2 has
alre ady shown that
the e dge -se t of D N , which is a
subse t of
the e dge -se t of D P , is
su#cie nt to
asse mble L k,j .
We now show that L k,j will
be absorbe d from C i by
one of i's
pare nts in D P
adde d to C i
due to pivoting,
irre spe ctive of
whe the r row k
be longe d to
the original Struct(L #,i ) or if it too was
adde d to C i
due to pivoting.
LU-pare nt(i) and
LU-pare nt(j).
We conside r two
case s: (1) g # h
and (2) g > h. If g # h, F g will
have both row k and column j and will absorb
the e le me nt
corre sponding to L k,j from C i . If g > h,
the n
the first condition for
the applicability of
The ore m 4.5 has
be e n
d. Now
we conside r two
furthe r
sce narios: (2a) At
le ast
one of i's
U-pare nts is in
the original Struct(L #,j ), or (2b)
none of i's
U-pare nts is in
the original Struct(L #,j ). In
case of (2a),
afte r pivoting,
at
le ast
one of i's
U-pare nts is in
the ne w Struct(U #,j ) and
the frontal matrix of this
U-pare nt will absorb column j from C i , including
the e ntry
corre sponding to L k,j . In
case (2b),
the se cond condition for
the applicability of
The ore m 4.5 has
be e n
d.
Finally,
whe the r row k was in
the original Struct(L #,i ) or was
adde d to C i
due to
the failure of a
U-de sce nde nt k, in its final location, k must
be gre ate r than h.
The re ason
is that if j # k # h
(i.e ., k's
ne w location is in
the e xte nde d
supe rnode h),
the n h
must
be an
LU-ance stor of i
be cause
implie s that
the re are both i L
#h and
iU #h in
the data-DAG. But that is not
possible be cause we are alre ady working
unde r
the assumption that
the LU-pare nt g of i is
gre ate r than h.
The re fore , k > h and
the third condition of
The ore m 4.5 has also
be e n
d. As a
re sult,
The ore m 4.5
would
have e nsure d that a
U-e dge iU #h is
pre se nt in D P to
asse mble column j from
Similarly,
we can
prove that
noe ntry
corre sponding to any U j,k will
be le ft
unasse mble d in C i .
4.6. Experimental results. In
se ctions 4.4 and 4.5,
we showe d how to supple
me nt
the e dge -se t of
the task-DAG to construct a data-DAG for
the unsym-
me tric multifrontal algorithm.

Table
showse xpe rime ntal
re sults of WSMP's
UNSYMMETRIC SPARSE MATRIX FACTORIZATION 545

Table
Time required for
and the number of edges in each DAG.
Matrix Symbolic Supplement-1 Supplement-2 Total Time
|ED P |
bbmat 1.7 6077 .06 6142 .00 6181 1.76 1.02
mil053 4.5 166154 .51 166154 .07 166154 5.08 1.00
mixtank 1.2 3203 .05 3203 .00 3203 1.25 1.00
pre2 6.4 317216 .84 320063 .16 320942 7.40 1.01
twotone .91 44856 .12 45866 .01 45918 1.04 1.02
wang3old .54 8450 .03 8450 .00 8450 0.57 1.00
imple me ntation of
the proce dure s to
ge ne rate the various DAGs.
Thre e DAGs
are conside re d in

Table

4.1:
the supe rnodal task-DAG T S ,
the supe rnodal data-DAG
D N for
unsymme tric multifrontal factorization without pivoting, and
the supe rnodal
data-DAG D P for
unsymme tric multifrontal factorization with pivoting.
The table shows
the time to
compute e ach of
the DAGs and
the numbe r
dge s in
the m for
the
matrice s in our
te st
suite .
T S is
compute d by
the basic symbolic factorization algorithm
de scribe d in
se ction
the re fore , t S is
the basic symbolic factorization
time .
We re fe r to
the proce ss of
computing D N from T S as
Supple me nt-1.
Supple me nt-1
che cks for
the first two
conditions of
The ore m 4.1 to find
the e dge s to
be adde d to E T S and
the n adds
outgoing
dge s from
supe rnode s without
LU-pare nts to
ld E D N .
Supple me nt-2
is
the proce ss that
addse dge s
base d on
the first two conditions of
The ore m 4.5 to
to
ld
The e xe cution
time of
Supple me nt-1 and
Supple me nt-2 is
de note d
by t 1 and t 2 ,
re spe ctive ly.
Note that not all
the e dge s in D N and D P
are ne ce ssary. For
the sake of com-
putational
spe e d,
Supple me nt-1 and
Supple me nt-2 in WSMP do not
che ck for all
the conditions of
The ore ms 4.1, 4.4, and 4.5
while addinge dge s.
The last conditions of all
thre e the ore ms
are skippe d.
Eve n if all conditions of
the se the ore ms
we re che cke d, not
all
the e dge s in
the re sulting data-DAGs may
be ne ce ssary.
The re fore , D N and D P
are not minimal data-DAGs for
unsymme tric multifrontal factorization.
Howe ve r, as

Table

4.1 shows,
the se DAGs do not
have many
more e dge s than T S for most
re al-life matrice s.
The ave rage fore xce sse dge s in
supe rnodal D P
ove r T S is only about 4%
for our
te st
suite .
We have shown that
the e dge s in
the task-DAGs
or T S
are
insu#cie nt to
dire ct
the data flow in
unsymme tric multifrontal factorization. On
the othe r hand,
the e dge s in D P
are su#cie nt,e ve n with pivoting.
The re fore ,
the numbe r
dge s in a truly minimal
supe rnodal data-DAG is
some whe re be twe e n
the numbe r
dge s in T S and in
the supe rnodal D P .
The e xpe rime ntal
re sults in

Table

4.1 show
that
the se two
rs
are usually fairly
close .
The table also shows that
the time re quire d to construct D N and D P is also small
compare d to
the basic symbolic
torization
time . Thus,
the me thodology
de scribe d in this
se ction for
the construction
of data-DAGs for
unsymme tric multifrontal factorization
nt in both
time and
the numbe r of
DAGe dge s. A comparison of
the

Table

4.1 with
WSMP's LU factorization
time give n in

Table

5.1 shows that
the total symbolic
time is usually significantly
le ss than
the nume rical factorization
time .
5. Implementation details of unsymmetric factorization. A
brie f
outline of
the unsymme tric multifrontal algorithm
base d on
the work of
Hadfie ld [20] and
Davis and Du# [9] is found in
se ction 4.2.
We now add
more de tails to it and
pre se nt
a
comple te algorithm that is
imple me nte d in WSMP. WSMP is
ge are d towards multiple
factorizations of
matrice s with
the same sparsity
patte rn but
di#e re nt
nonze ro
value s.
The re fore , symbolic
phase is
pe rforme d only
once and its output is
use d in
all
subse que nt
nume rical
ve n with
di#e re nt pivot
se que nce s
re sulting
from
di#e re nt
nume rical
value s.
A
fundame ntal
data-structure in our
unsymme tric multifrontal algorithm is
the frontal matrix. A frontal matrix is
associate d
withe ach
supe rnode .

Figure

5.1 shows
the organization of a typical frontal matrix for a
supe rnode
The core of this frontal matrix is
a| Struct(L #,q
)| -| Struct(U q,#
)| portion,
re
are pre dicte d by
the symbolic factorization. In
the abse nce of pivoting,
the first r - q columns of this matrix would
be factore d and would
be save d as parts of U and L,
re spe ctive ly.
The re maining trailing
submatrix would
constitute the contribution matrix
whose conte nts would
be absorbe d
into
the frontal
matrice s of
the pare nts of g in D P .
Extra rows
F [q:r]
s [q:r]
Extra columns
r
pivot block
pivot block
Fig . 5.1. Organization of a typical frontal matrix for a supernode r]). The p failed
pivots from the LU-children of the supernode are appended at the beginning of the frontal matrix
and the extra
row s and columns inherited from U- and L-descendents, respectively, are appended at
the end.
UNSYMMETRIC SPARSE MATRIX FACTORIZATION 547
In
the pre se nce of
nume rical
pivoting,e xtra pivots as
we ll as
othe r rows and
columns may
be adde d to
the frontal matrix
de pe nding on
the labe ls and pivot
failure s
of
the childre n of g in D P . Extra pivots (row-column pairs with
the same indice s)
are adde d to F g if
some of
the pivots of g's
LU-childre n fail to satisfy
the pivoting
crite rion.
The LU-childre n of g
the mse lve s may
have inhe rite d
some or all of
the se faile d pivots from
the ir own
LU-childre n.
The re fore ,
faile d pivots from any of
the LU-de sce nde nts of g
cane nd up in its frontal matrix. If p such pivots
are adde d,
the n
the size of
the pivot block
incre ase s from r - q
The frontal matrix F g can similarly
inhe rite xtra rows
corre sponding to
faile d
pivots in its
U-de sce nde nts
whose LU-pare nts
are gre ate r than g
ande xtra columns
corre sponding to
faile d pivots in its
L-de sce nde nts
whose LU-pare nts
are gre ate r than
g.
Irre spe ctive of
the ir
ne w
indice s,
the se e xtra rows and columns
are always ap-
pe nde d at
the e nd of
the original rows and columns of F g and a
d list of
the ir
indice s is
maintaine d
ate ach
supe rnode .
Eve ntually,
the se are asse mble d into
the e xtra pivots of
the frontal
matrice s of
the LU-pare nts of
the supe rnode s
whe re the se pivots
faile d.
The row and column
structure s
pre dicte d by symbolic factorization
are intact for
future factorizations of
matrice s with
the same nonze ro
patte rn.
The additions to
the se structure s
due to pivoting, which
de pe nd on
the nonze ro
value s in
the matrix
be ing
factore d,
are maintaine d
se parate ly and
are discarde d
be fore e ach
ne w factorization.
The availability of a static data-DAG D P that is
su#cie nt for handling an arbitrary
amount of dynamic pivoting is critical to our
imple me ntation of
the unsymme tric
multifrontal algorithm.

Figure
give s a
high-le ve l
pse udocode of our factorization
algorithm.
The algorithm starts with
the root
supe rnode of task- and data-DAGs. At
any
supe rnode , first, it
re cursive ly factors all
the unfactore d
childre n of that
supe r-
node .
The n it looks at
the faile d pivots (if any) of its
childre n to
figure out
the numbe r and
indice s of
the e xtra rows, columns, and pivots, if any, and accordingly
allocate s a frontal matrix of
the appropriate size . In
the ne xt
ste p,
the contribution
from
the original
coe #cie nt matrix and
the contribution
matrice s of
the curre nt
supe rnode 's
childre n
are accumulate d in
the appropriate locations
inside the frontal
matrix. Finally,
the algorithm
proce e ds to factor
the pivot block of
the frontal matrix
and
update s
the re mainde r of
the frontal matrix.
The le ading
succe ssfully
factore d
rows and columns
are save d as portions of U and L for
use during triangular
solve s.
The re maining contribution matrix
ve ntually
asse mble d into
the frontal
matrice s
of its
pare nts and is
re le ase d by
the last
pare nt to pick up its contribution.
The frontal matrix of
the LU-pare nt of a
supe rnode picks up all its
faile d pivot
row-column pairs as
we ll as
the e ntire trailing submatrix of its contribution matrix
with row and column
indice s
gre ate r than
ore qual to
the first
inde x of
the pare nt
supe rnode .
The re maining rows and columns of a
supe rnode 's contribution matrix
are asse mble d into
the frontal
matrice s of its L- and
U-pare nts in D P . It is
possible for
more than
one L- or
U-pare nts' frontal
matrice s to
have the same row or column
indice s in common with
the child's contribution matrix.
Howe ve r,e ache le me nt of
a contribution matrix must
be adde d
intoe xactly
one frontal matrix.
Some simple bookke e ping to
track of rows and columns that
have be e n
asse mble d
su#ce s to
e nsure this condition for
the re lative ly
fe w rows and columns that
have the pote ntial
to
be copie d into
the frontal
matrice s of
multiple L- and
U-pare nts,
re spe ctive ly.

Figure

5.2 and
the de scription in this
se ction show that WSMP's
unsymme tric
multifrontal algorithm is fairly straightforward to
imple me nt.
The static task- and
data-DAGs
compute d during
the symbolic
phase and
the use of
re cursion
make the
548 ANSHUL
function uns mf (root) {
/* 1.
Re cursive calls to root's
childre n */
fore ach child k of root in T S do
if not
alre ady
proce sse d k then
Call uns mf (k);
Flag
supe rnode k as
alre ady
proce sse d;
end for
/* 2.
Colle ct pivoting info to
de te rmine size of F root */
fore ach child k of root in D P do
if k is an L-child then
if k has
faile d pivots then
Add
the m to
the sorte d list of F root
'se xtra columns;
hase xtra columns then
Add
those whose LU-pare nt is
gre ate r than root to
the sorte d list of F root
'se xtra columns
while che cking for
duplicate s;
else if k is a U-child then
if k has
faile d pivots then
Add
the m to
the sorte d list of F root
'se xtra rows;
hase xtra rows then
Add
those whose LU-pare nt is
gre ate r than root to
the sorte d list of F root
'se xtra rows
while che cking for
duplicate s;
else if k is an LU-child then
if k has
faile d pivots then
Add
the m to
the sorte d list of F root
'se xtra pivots;
hase xtra columns then
Add
those whose LU-pare nt is
gre ate r than root to
the sorte d list of F root
'se xtra columns
while che cking for
duplicate s;
hase xtra rows then
Add
those whose LU-pare nt is
gre ate r than root to
the sorte d list of F root
'se xtra rows
while che cking for
duplicate s;
end for
/* 3.
Initialize root's frontal matrix */
Allocate F root of
appropriate size and fill it with
ze ros;
Populate F root
withe ntrie s from A
corre sponding to
supe rnode root;
/* 4.
Asse mbly from
childre n's contribution
matrice s into F root */
fore ach child k of root in D P do
Copy
appropriate contribution from C k into F root ;
if root is
the last
pare nt of k to pick up C k 's contribution then
Fre e the space occupie d by C
end for
5.
Nume rical factorization */
Factor
the pivot block of F root and
update the trailing part to
ld C root ;
function uns mf.
Fig . 5.2. A simple and e#cient unsymmetric multifrontal algorithm.
UNSYMMETRIC SPARSE MATRIX FACTORIZATION 549
nume rical factorization algorithm much
simple r to
de scribe and
imple me nt than
the e arlie r
de scriptions of
the unsymme tric
patte rn multifrontal algorithm in [20] and [9].
Othe r than UMFPACK [8], WSMP is
the only
sparse solve r
available that is
base d on
an
unsymme tric
patte rn multifrontal algorithm. It is also
the first such
paralle l
solve r
available for
ge ne ral
use . Although
Hadfie ld [20]
provide de xpe rime ntal
re sults from a
paralle l
imple me ntation on
the nCUBE, a practical
paralle l
solve r did not
re sult from
thate #ort.
The algorithm of

Figure

5.2 is not only
re lative ly
simple in
de scription but is
also computationally
le an
be cause it
the none sse ntial non-floating-point
ope rations and can
handle pivot
failure s
fairlye #cie ntly. It is also
note worthy that
for structurally
symme tric
matrice s,
the algorithm in

Figure

5.2 naturally
re duce s
to a
symme tric-patte rn multifrontal algorithm
guide d by
ane limination
tre e , which
re place s both T S and D P .
Othe r than a
state me nts
fore ach
supe rnode ,
the re is no
ove rhe ad in using this algorithm for structurally
symme tric
matrice s.
5.1. Experimental results.
We now
compare the unsymme tric LU factoriza-
tion
time of WSMP with that of
thre e state -of-the -art multifrontal
sparse solve rs,
name ly, MUMPS
ve rsion 4.1.6 [4, 5], MA41 [2, 3], and UMFPACK
ve rsion 3.2 [8]. A
de taile d
comparative study that
include s
more solve rs can
be found in [18].
The software
compare d in this
se ctione mploy
di#e re nt variants of
the multifrontal
me thod.
MUMPS contains a
symme tric-patte rn multifrontal factorization
code base d on
the classical multifrontal algorithm [14]. MA41, in
some se nse , is a hybrid
be twe e n sym-
me tric and
unsymme tric
patte rn multifrontal
solve rs. It
use s
ane limination
tre e to
guide factorization, but
the frontal
matrice s
are prune d of
all-ze ro rows and columns.
UMFPACK 3.2 contains a variation of
the unsymme tric-patte rn multifrontal algorithm
[9] that
use s
ane limination
tre e de rive d from
the structure of A # A.
Apart from
the factorization algorithm,
the re are othe r significant
di#e re nce s
among
the four
software package s that
a#e ct
the ir
pe rformance . First,
the y
use diffe
re nt
sche me s for
fill-re ducing
ring. By
de fault, WSMP
use s a
symme tric
pe rmu-
tation
base d on a
ne ste d-disse ction
ring [17]
compute d on
the structure of A+A # .
MUMPS and MA41
use a
symme tric
pe rmutation
base d on
the approximate mini-
mum
de gre e (AMD) algorithm [1]
applie d to
the structure of A+A # . UMFPACK
use s
a column AMD algorithm [10] to
pre pe rmute only
the columns of A and
compute s a
row
pe rmutation
base d on
nume rical and sparsity
crite ria during factorization.
The se cond
di#e re nce is
the use of a maximal matching algorithm [13] to
pe rmute the rows
of
the coe #cie nt matrix to
maximize the product of
the magnitude s of its diagonal
e ntrie s. As shown in [6, 18], this can
a#e ct factorization
time s
be cause it
change s
the amount of structural
symme try and
the amount of
nume rical pivoting during factor-
ization. WSMP
use s this
pre proce ssing on all
matrice s, MUMPS and MA41
use it
only if
the structural
symme try in
the original matrix is
le ss than 50%, and UMF-
doe s not
use it at all.
The third
di#e re nce is that WSMP
re duce s
the coe #cie nt
matrix into a block-triangular form,
while MUMPS, MA41, and UMFPACK 3.2 do
not.

Table

5.1 shows
nume rical factorization
time s and
ope ration counts of MUMPS,
MA41, UMFPACK, and WSMP run with
the options in MUMPS, MA41, and WSMP
change d to
minimize the di#e re nce s
be twe e n
the code s
othe r than
the factorization
algorithm.
We switche d o#
the pe rmutation to a
he avy-diagonal form and
the associate
d scaling in MUMPS, MA41, and WSMP. For WSMP,
inste ad of its
de fault
ne ste d-disse ction
ring,
we use d an
approximate minimum fill
orde ring, which is
ve ry similar to AMD.
Eve n with
the se change s,
di#e re nce s
re main
be twe e n
the four

Table
LU factorization times and operation counts of MUMPS, MA41, UMFPACK 3.2, and WSMP
w ith similar permutation options. The best time is in boldface and the second best time is underlined.
3.2W SMP
Matrix time ops time ops time ops time ops
af23560 3.89 2.56 3.58 2.54 8.59 3.46 4.06 3.22
bbmat 54.3 41.6 56.3 41.1 78.7 39.1 27.6 21.5
e40r0000 4.93 2.53 3.63 1.58 6.23 2.17 0.80 .419
ecl32 64.2 64.6 67.1 64.4 191. 112. 139. 77.6
fidap011 8.58 7.01 8.79 6.96 17.0 8.51 6.51 5.74
mil053 43.5 31.8 40.0 31.8 107. 46.2 28.2 20.8
mixtank 151. 141. 152. 64.1 363. 243. 76.3 64.6
nasasrb 12.8 9.45 11.9 9.43 55.9 28.2 10.4 8.78
pre2 fail fail fail fail fail fail 346. 301.
raefsky3 4.44 2.90 3.88 2.90 16.0 7.87 4.88 4.17
raefsky4 107. 74.4 92.9 44.7 25.0 12.9 43.4 22.5
twotone 56.5 38.3 37.6 31.8 30.1 10.8 2.87 1.49
wang3old 72.9 57.8 57.7 51.0 40.6 24.2 45.8 32.3
wang4 11.8 10.5 12.2 10.5 53.4 30.7 8.84 7.94
code s. For
instance , MUMPS, MA41, and WSMP first
pe rmute the matrix such
that it has a diagonal of structural
nonze ros. This initial
pe rmutation is
the same for
MUMPS and MA41
be cause both
use the same code to
compute it.
Howe ve r, it can
be di#e re nt for WSMP.
The pivoting
strate gy of UMFPACK
base d on row
is
inhe re ntly
di#e re nt from
the symme tric
inte rsupe rnode pivoting
strate gy
use d in
MUMPS, MA41, and WSMP. WSMP's algorithms work only with a
pe rmutation to
the block-triangular form, which is not
imple me nte d in MUMPS, MA41, and UMF-
PACK.
Howe ve r, with
the e xce ption of comp2c,
the e #e ct of block-triangularization
on
the ope ration count for factorization is insignificant, if any. As a
re sult of
the se di#e re nce s and
due to
the fact that MUMPS may
pe rform
more ope rations than
ne c-
e ssary on structurally
unsymme tric
matrice s,
the factorization
ope ration counts for
the four
code s in

Table
are di#e re nte ve n with a similar
ring algorithm for
fill-re duction.
In

Table

5.1,
the faste st factorization
time fore ach matrix is in
boldface and
the se cond
faste st
time is
unde rline d. Although
di#e re nce s
othe r than
the factorization
algorithm
itse lf
a#e ct
the pe rformance of
the se code s, it
ise asy to
se e the broad
picture thate me rge s from

Table

5.1. Most of
the boldface e ntrie s
are in
the WSMP column
and most of
the unde rline de ntrie s
are in
the MA41 column. For many
matrice s,
the e #e ct of
the algorithmic
choice s of
the software ise vide nt in
the factorization
statistics in

Table

5.1. MUMPS usually
re quire s
more floating-point
ope rations for
factorization than MA41 and WSMP
be cause it
use s artificially
symme trize d frontal
matrice s
padde d with
ze ros. For
the same re ason, UMFPACK is
faste r than MUMPS
UNSYMMETRIC SPARSE MATRIX FACTORIZATION 551
for
ve ry
unsymme tric
matrice s (such as
baye r01,
one tone 2, and
twotone
howe ve r,
it is
slowe r for
matrice s with
more structural
symme try (such as fidap011, mil053,
and wang4) partly
be cause it
use s a
fill-re ducing
pe rmutation on
the columns of
the coe #cie nt matrix
be fore starting LU factorization. MA41
o#e rs a significant
improve me nt
ove r MUMPS for
matrice s with a
ve ry
unsymme tric
structure , such as
comp2c,
one tone 1, and
twotone . It
se e ms that MA41's
me chanism for finding all-
ze ro rows and columns incurs a slight
ove rhe ad that it cannot
o#se t for
matrice s
with a
ne arly
symme tric
structure (such
ase cl32, fidap011, and wang4), for which
it is
some what
slowe r than MUMPS. It is
cle ar from

Table

5.1 that WSMP has
the smalle st
ove rall factorization
time se ve n
its
de fault options
are modifie d to
compare it with
the othe r
solve rs. With its
de fault options, WSMP's factorization
time s
are usually much
smalle r [18] than
those shown in

Table

5.1.
6. Concluding remarks. This
pape r
de scribe s
sparse unsymme tric symbolic
and
nume rical factorization algorithms that
improve pre vious similar algorithms. Our
phase , in particular, is
more powe rful than
othe rs
de scribe d in
the lite rature . It
ine xpe nsive ly
compute s
minimale limination
structure s that
are transitive re ductions of
the uppe r and
lowe r triangular factors of
the original
coe f-
ficie nt matrix. In addition, it
compute s
ne ar-minimal
data-de pe nde ncy DAGs for
unsymme tric multifrontal factorization with and without pivoting. A data-DAG that
has only a slightly
highe r
ofe dge s than a minimal task-DAG and that is
capable ofe xpre ssing all
possible data-de pe nde ncie s in
the face of dynamic pivoting
is a
fe ature of our symbolic
phase .
We show how this data-DAG aids in a
high-pe rformance imple me ntation of
the unsymme tric multifrontal LU factorization
algorithm. This factorization algorithm is not only
faste r than
othe r
sparse LU factorization
algorithms but is also
simple r than
the unsymme tric multifrontal algorithms
de scribe d
pre viously in
the lite rature .
Our algorithms do not
introduce additional
ove rhe ads
while factoring
matrice s
with a
symme tric
nonze ro
patte rn.
Whe n
pre se nte d with a
sparse matrix with a
symme tric
structure , both
the symbolic and
the nume rical factorization algorithms
and
the data-structure s
ge ne rate d by
the m
grace fully transform into
the ir
symme tric
counte rparts without
re quiring any significant amount
ofe xtra
proce ssing or
storage .
In a
distribute d-me mory
paralle l
imple me ntation of
unsymme tric
sparse LU
the e dge s of
the data-DAG
conne cting tasks
mappe d onto
di#e re nt
proce sse s
de te rmine the inte rproce ss communication
patte rn.
The static and
ne ar-
minimal
nature of
the data-DAG
use d in our algorithms would
be e xtre me ly
use ful for
pote ntial
paralle l
imple me ntations of
unsymme tric multifrontal factorization,
re changing
the data-DAG dynamically could
be cumbe rsome and
ine #cie nt and
the unne ce ssary
dge s could
incre ase synchronization and communication
ove rhe ads.

Acknowledgments

.
The author
wishe s to thank
Andre w Conn,
Fre d Gustavson

Jose ph Liu, Sivan
Tole do, and
the anonymous
re fe re e s
fore xtre me ly
use ful
comme nts
one arlie r drafts of this
pape r.



--R

An approximate minimum degree ordering algorithm
Vectorization of a multiprocessor multifrontal code
Memory management issues in sparse multifrontal methods on multiprocessors
A fully asynchronous multifrontal solver using distributed dynamic scheduling
Multifrontal parallel distributed symmetric and unsymmetric solvers
Analysis and comparison of tw sparse solvers for distributed memory computers
The influence of relaxed supernode partitions on the multifrontal method

An unsymmetric-pattern multifrontal method for sparse LU factorization
A Column Approximate Minimum Degree Ordering Algorithm
A supernodal approach to sparse partial pivoting
Direct Methods for Sparse Matrices
On algorithms for permuting large entries to the diagonal of a sparse matrix
The multifrontal solution of unsymmetric sets of linear equations
Exploiting structural symmetry in unsymmetric sparse symbolic factorization
Elimination structures for unsymmetric sparse LU factors
Fast and e
Recent advances in direct methods for solving unsymmetric sparse systems of linear equations
Highly scalable parallel algorithms for sparse matrix factorization
On the LU Factorization of Sequences of Identically Structured Sparse Ma- tricesw ithin a Distributed Memory Environment
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
--TR

--CTR
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 matrix factorization;parallel sparse solvers;multifrontal methods