Systems Optimization Laboratory
Stanford, CA 94305-4121 USA
LUSOL: Sparse LU for Ax = b
- CONTRIBUTORS: Philip Gill, Walter Murray, Margaret Wright,
Michael O'Sullivan, Kjell Eikland, Yin Zhang, Nick Henderson
- CONTENTS: A sparse LU factorization for square and rectangular
matrices A, with Bartels-Golub-Reid updates for column replacement
and other rank-1 modifications.
Typically used for a sequence of linear equations as in the
Solve Ax = b and/or A'y = c
Replace a column of A
Repeat with different b, c
The matrix A may have any shape and rank.
Rectangular LU factors may be used to form a sparse null-space matrix operator.
Special feature 1:
Three sparse pivoting options in the Factor routine:
Threshold partial pivoting (TPP)
Threshold rook pivoting (TRP)
Threshold complete pivoting (TCP)
All options choose row and column permutations as they go,
balancing sparsity and stability according to different rules.
TPP is normally most efficient for solving Ax = b.
TRP and TCP are rank-revealing factorizations.
In practice, TRP is an effective method for estimating rank(A).
TCP tends to be too dense and expensive to be useful.
Special feature 2:
Multiple update routines:
Add, delete, or replace a column of A
Add, delete, or replace a row of A
Add a general (sparse) rank-1 matrix to A
The condition of L is controlled throughout, so that U
tends to reflect the condition of A.
This is essential for subsequent Bartels-Golub-type updates
(which are implemented in a manner similar to John Reid's LA05
and LA15 packages in the HSL library).
If a fresh factorization is thought of as A = LDU
(with unit diagonals on L and U), then TRP and TCP
control the condition of both L and U. This is why
they have rank-revealing properties.
LUSOL is the basis factorization package (BFP) for
Factor: No special handling of dense columns.
Solve: No special treatment of sparse right-hand sides.
Documentation: No user's manual. Primary documentation is
in-line comments within the f77 source code,
and README files within the src, csrc, test, matlab, cmex directories.
J. K. Reid (1982).
A sparsity-exploiting variant of the Bartels-Golub
decomposition for linear programming bases,
Mathematical Programming 24, 55-69.
P. E. Gill, W. Murray, M. A. Saunders and M. H. Wright (1987).
Maintaining LU factors of a general sparse matrix,
Linear Algebra and its Applications 88/89, 239-270.
P. E. Gill, W. Murray and M. A. Saunders (2005).
SNOPT: An SQP algorithm for large-scale constrained optimization,
SIGEST article, SIAM Review 47(1), 99-131.
(See sections 4 and 5.)
- 01 Feb 2008:
First downloadable version, including f77 source code, C translation
by Kjell Eikland, and two cmex implementations by Mike O'Sullivan
and Yin Zhang respectively.
- 06 Dec 2009:
Example Matlab files for forming a well-conditioned nullspace operator Z
from LUSOL's LU factors of a sparse rectangular matrix, and applying
it to a given vector or matrix: Z*v or Z*V.
Based on 01 Feb 2008 lusol.zip.
- 17 Jan 2011:
New implementation of LUSOL mex interface and M-files by
Nick Henderson, ICME, Stanford.