Algorithm 8xx: UMFPACK V4.1, an
unsymmetricpattern multifrontal method
Timothy A. Davis*
May 6, 2003
Abstract
An ANSI C code for sparse LU factorization is presented that combines a column
preordering strategy with a rightlooking unsymmetricpattern multifrontal numerical
factorization. The preordering and symbolic analysis phase computes an upper bound
on fillin, work, and memory usage during the subsequent numerical factorization.
Usercallable routines are provided for ordering and analyzing a sparse matrix, com
puting the numerical factorization, solving a system with the LU factors, transposing
and permuting a sparse matrix, and converting between sparse matrix representations.
The simple user interface shields the user from the details of the complex sparse fac
torization data structures by returning simple handles to opaque objects. Additional
usercallable routines are provided for printing and extracting the contents of these
opaque objects. An even simpler way to use the package is through its MATLAB in
terface. [ FPPACK is incorporated as a builtin operator in MATLAB 6.5 as x = A\b
when A is sparse and unsymmetric.
Categories and Subject Descriptors: G.1.3 \ii'i idil AnIi1.i.]: niin.!i iil Linear Al
gebra linear systems (direct methods), sparse and very large systems G.4 [i'I, !i.!e ii, of
C'! 11l,,1 iil]: Mathematical Software , d',i .:I/,, analysis, ti .:' '. ..
General terms: Algorithms, Experimentation, Performance.
Keywords: sparse !iii. !.: !ii!r1 lie matrices, linear equations, multifrontal method, order
ing methods.
1 Overview
UMIFPACK Version 4.1 is a code written in ANSI C for the direct solution of systems of
linear equations, Ax = b, where A is sparse and tii.: iiii l!i The matrix PAQ, PRAQ,
or PR AQ is factorized into the product LU. The column ordering Q is selected to give a
*Dept. of Computer and Information Science and En11li1.. h 1i, Univ. of Florida, Gainesville, FL, USA.
email: davis@cise.ufl.edu. http://www.cise.ufl.edu/~davis. This work was supported by the National Science
Foundation, under grants ASC9111263, DMS9223088, and DMS(i2'i..12'7 Portions of the work were done
while on sabbatical at Stanford University and Lawrence Berkeley National Laboratory (with funding from
Stanford University and the SciDAC program).
good a priori upper bound on fillin and then refined during numerical factorization (while
preserving the upper bound on fillin). The row ordering P is selected during numerical
factorization to maintain numerical stability and to preserve sparsity. The diagonal matrix
R scales the rows of the matrix A. The method is a combination of a column preordering
and symbolic ,i!ll.i., phase based on COLAII) or AIlI) and a numerical factorization
phase based on a modification of the uiii. inliii1 i:pattern multifrontal method (UMFPACK
Version 2, or 1. 38, which does not have a preordering and symbolic iiii1,.i., phase). The
methods used by UMFPACK V4.1 are discussed in a companion paper [1]. The sparse matrix
A factorized by UMIFPACK can be real or complex, square or rectangular, and singular or
nonsingular (or any combination).
U;I FPACK ,in 1 .1. the matrix and automatically selects between three ordering strate
gies, described below.
unsymmetric: A column preordering is computed by a modified version of COLA ,II).
The method finds a permutation Q that limits the fillin for any subsequent choice of
P via partial pivoting. This modified version of COLAMII) also computes the column
elimination tree and a depthfirst postordering of the tree. During factorization, the
column preordering can be modified. Columns within a single supercolumn can be
reshuffled, to reduce fillin. Threshold partial pivoting is used with no preference given
to the diagonal entry. Within a given pivot column j, an entry aij can be chosen if
aij > 0.1 max a,j Among those numerically acceptable entries, the sparsest row i is
chosen as the pivot row.
symmetric: The column ordering is computed from AMII) applied to the pattern of
A+AT, followed by a postordering of the supernodal elimination tree. No modification
of the column preordering is made during numerical factorization. Threshold partial
pivoting is used, with a strong preference given to the diagonal entry. The diagonal
entry is chosen if ajj > 0.001 max aj Otherwise, a sparse row is selected, using the
same method used by the tii!.: iiiiil1 ii: strategy (with a relative threshold of 0.1).
2by2: A row permutation P2 is found which attempts to reduce the number of small
diagonal entries of P2A. If aii is numerically small, the method attempts to swap
two rows i and j, such that both aij and aji are large. Once these rows are swapped
they remain in place. This does not guarantee a zerofree diagonal, but it does tend
to preserve the symmetry of the nonzero pattern. \_':i, the symmetric strategy (see
above) is applied to the matrix P2A.
2 MATLAB Interface
The MIATLAB interface to UMFPACK provides a replacement for ,IATLAB's LU routine,
and the forward slash and backlash matrix operators. It is typically much faster than the
builtin routines in MAI TLAB Version 6.0, uses less memory, and returns sparser LU factors.
II ATLAB 6.5 includes UMIFPACK Version 4.0 as a builtin routine. The interface provided
here also allows access to U,IFPACK's preordering and symbolic !il1:.i., phase.
3 ANSI C Interface
The ANSI C U, IFPACK library consists of 31 usercallable routines and one include file.
Twc!l .i:, !i of the routines come in four versions: real or complex (both double precision),
and int or long integers. Only the double / int version is described here; the other versions
are analogous. UMIFPACK requires the BLAS (which perform dense matrix operations) for
best performance, but can be used without the BLAS. Fi e primary UMIFPACK routines
are required to solve Ax = b:
umfpack_di_symbolic: Preorders the columns of A, finds the supernodal column elim
ination tree, and postorders the tree. Returns an opaque Symbolic object as a void *
pointer that can be used by umfpacknumeric to factorize A or any other matrix with
the same nonzero pattern as A. Computes upper bounds on the nonzeros in L and U,
the floatingpoint operations required, and the memory usage of umfpack_di_numeric.
umfpack_dinumeric: \niiri i, illy factorizes a sparse matrix PAQ, PRAQ, or PR AQ
into the product LU, using the Symbolic object. Returns an opaque Numeric object
as a void pointer.
umfpack_disolve: Solves a sparse linear system (Ax = b, ATx = b, or systems
involving just L or U), using the Numeric object. Performs iterative refinement with
sparse backward error.
umfpack_free_symbolic: Frees the Symbolic object.
umfpack_freenumeric: Frees the Numeric object.
The matrix A is represented in compressed column form:
int Ap [n+1] ;
int Ai [nz] ;
double Ax [nz]
The row indices and numerical values of entries in column j are stored in Ai[Ap[j] ...
Ap[j+1]1] and Ax[Ap[j] ... Ap[j+1]11, respectively. This simple program illustrates the
basic usage of UMIFPACK:
#include
#include "umfpack.h"
int n = 5 ;
int Ap [ ] = {0, 2, 5, 9, 10, 12} ;
int Ai [ ] = { 0, 1, 0, 2, 4, 1, 2, 3, 4, 2, 1, 4 ;
double Ax [ ] = {2., 3., 3., 1., 4., 4., 3., 1., 2., 2., 6., 1. ;
double b [ ] = {8., 45., 3., 3., 19. ;
double x [5] ;
int main (void)
{
double *null = (double *) NULL ;
int i ;
void *Symbolic, *Numeric ;
(void) umfpack_di_eimbJclic (n, n, Ap, Ai, Ax, rSimbJclic, null, null)
(void) umfpack_di_numeric (Ap, Ai, Ax, Symbolic, &Numeric, null, null)
umfpack_di_free_symbolic (&Symbolic) ;
(void) umfpack_di_solve (UMFPACK_A, Ap, Ai, Ax, x, b, Numeric, null, null)
umfpack_di_free_numeric (&Numeric) ;
for (i = 0 ; i < n ; i++) printf ("x [Ed] = %g\n", i, x [i]) ;
return (0) ;
The Ap, Ai, and Ax I! i,., represent the matrix
2 3 0 0 0
3 0 406
A= 0 1 3 2 0
0 0 100
0 4 201
and the solution is x = [12345]T. Additional routines are provided for:
Changing default parameter settings, and for providing a different column preordering.
Converting a matrix in triplet form to compressed column form, and visa versa. The
triplet form is a simpler data structure for the user to manipulate. It consists of three
~!Ii,,.I that hold the row index, column index, and numerical value of each entry in
matrix.
Transposing and optionally permuting a compressed column form matrix.
Getting the contents of the opaque Symbolic and Numeric objects, saving them to a
file, and loading them from a file.
Printing and verifying control parameters, statistics, sparse matrices, Symbolic and
Numeric objects, permutation vectors, and dense vectors.
In addition to appearing as a Collected Algorithm of the ACM, UIFPACK is available
at http://www.cise.ufl.edu/research/sparse. The package includes a user guide that provides
full details on how to install the package, how to use the ATLAB interface, and how to
use the ANSI C interface. A basic Fortran interface is also provided.
References
[1] T. A. Davis. A column preordering strategy for the uii:. iiiiirl! iipattern multifrontal
method. ACM Trans. Math. Sof,.. 2003 (under submission). Also TR03006 at
www.cise.ufl.edu/techreports.
