Systems Optimization Laboratory
Stanford, CA 943054121 USA

LUSOL: Sparse LU for Ax = b
 AUTHOR:
Michael Saunders
 CONTRIBUTORS: Philip Gill, Walter Murray, Margaret Wright,
Michael O'Sullivan, Kjell Eikland, Yin Zhang, Nick Henderson, Ding Ma
 CONTENTS: A sparse LU factorization for square and rectangular
matrices A, with BartelsGolubReid updates for column replacement
and other rank1 modifications.
Typically used for a sequence of linear equations as in the
simplex method:
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 nullspace 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 rankrevealing factorizations.
In practice, TRP is an effective method for estimating rank(A).
TCP tends to be too dense and expensive to be useful, although
MINOS and SNOPT switch from TPP to TRP and even TCP if necessary
in case of persistent numerical difficulty.
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) rank1 matrix to A
Numerical stability:
LUSOL maintains LU factors with row and column permutations P, Q
such that A = LU with PLP' lower triangular (with unit diagonals)
and PUQ upper triangular.
The condition of L is controlled throughout by maintaining
Lij <= factol (= 10 or 5 or 2 or 1.1, ...), so that U
tends to reflect the condition of A.
This is essential for subsequent BartelsGolubtype 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 PLP' and PUQ), then TRP and TCP
control the condition of both L and U by maintaining
Lij <= factol and Uij <= factol, so that D
reflects the condition of A. This is why TRP and TCP have
rankrevealing properties.
Proven applications:
LUSOL is the basis factorization package (BFP) for
MINOS,
SQOPT,
SNOPT,
lp_solve,
AMPL/PATH,
GAMS/PATH.
Shortcomings:
Factor: No special handling of dense columns.
Solve: No special treatment of sparse righthand sides.
Documentation: No user's manual. Primary documentation is
inline comments within the f77 source code (and the more recent
f90 version).
 REFERENCES:
J. K. Reid (1982).
A sparsityexploiting variant of the BartelsGolub
decomposition for linear programming bases,
Mathematical Programming 24, 5569.
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, 239270.
P. E. Gill, W. Murray and M. A. Saunders (2005).
SNOPT: An SQP algorithm for largescale constrained optimization,
SIGEST article, SIAM Review 47(1), 99131.
(See sections 4 and 5.)
DOWNLOADS:
 01 Feb 2008:
lusol.zip
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:
nullspaceLUSOL.zip
Example Matlab files for forming a wellconditioned 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:
https://github.com/nwh/lusol_mex
New implementation of LUSOL mex interface and Mfiles by
Nick Henderson, ICME, Stanford.
 31 Jul 2013:
https://github.com/nwh/lusol
Second new implementation of LUSOL interface and Mfiles by
Nick Henderson, ICME, Stanford. This is based on a new f90 version of
the main factorization and columnreplace routines (lu1fac, lu6sol, lu8rpc).
In lu1fac the efficiency of Threshold Rook Pivoting is significantly
improved (Ding Ma and Michael Saunders, MSandE and ICME, Stanford).
Note: The github distribution includes 64bit binaries for Linux and Mac OS X.
