Development and history of sparse direct methods

Iain S. Duff
iain.duff@stfc.ac.uk

STFC Rutherford Appleton Laboratory

Oxfordshire, UK.
and
CERFACS, Toulouse, France
Homepage: http://www.numerical.rl.ac.uk/people/isd/isd.html

The SIAM Conference on Applied Linear Algebra. October 26-29, 2009. Monterey, California. – p.1/34
 Outline

Ancient times

The SIAM Conference on Applied Linear Algebra. October 26-29, 2009. Monterey, California. – p.2/34
 Outline

Ancient times
The birth

The SIAM Conference on Applied Linear Algebra. October 26-29, 2009. Monterey, California. – p.2/34
 Outline

Ancient times
The birth
Infancy

The SIAM Conference on Applied Linear Algebra. October 26-29, 2009. Monterey, California. – p.2/34
 Outline

Ancient times
The birth
Infancy
Turbulent teens

The SIAM Conference on Applied Linear Algebra. October 26-29, 2009. Monterey, California. – p.2/34
 Outline

Ancient times
The birth
Infancy
Turbulent teens
Maturity

The SIAM Conference on Applied Linear Algebra. October 26-29, 2009. Monterey, California. – p.2/34
 Outline

Ancient times
The birth
Infancy
Turbulent teens
Maturity
Senility?

The SIAM Conference on Applied Linear Algebra. October 26-29, 2009. Monterey, California. – p.2/34
 Direct methods
Gaussian Elimination
PAQ → LU
Permutations P and Q chosen to preserve sparsity and maintain stability
L : Lower triangular (sparse)
U : Upper triangular (sparse)

SOLVE:
Ax = b
by
Ly = Pb
then
UQT x = y

The SIAM Conference on Applied Linear Algebra. October 26-29, 2009. Monterey, California. – p.3/34
 Ancient times
Large systems were solved in ancient times by iterative methods,
principally relaxation methods (Richardson, Southwell, Fox )
and using SOR and related techniques ... thesis of David Young (1950),
book by Varga (1962).
The earliest references from standard sparse matrix books date to the early
1950s but these papers are principally concerned with graph theory and
combinatorics and don’t really reference or use sparse matrices or even
sparse data structures (König, Woodbury, Hall).

The SIAM Conference on Applied Linear Algebra. October 26-29, 2009. Monterey, California. – p.4/34
Ancient times

The SIAM Conference on Applied Linear Algebra. October 26-29, 2009. Monterey, California. – p.5/34
 Ancient times
One of the first papers to directly link graphs with sparse elimination was
due to Seymour Parter in 1961.
Major application areas were the simplex method and backward difference
formulae for solving stiff ODEs.
Other areas strongly represented in ancient times were:
Power systems (Bill Tinney)
Electrical networks and circuit design (Bob Brayton, Gary Hachtel,
Gabriel Kron, Alberto Sangiovanni-Vincintelli, Donald Steward)
Chemical engineering (Roger Sargent and Art Westerberg)

The SIAM Conference on Applied Linear Algebra. October 26-29, 2009. Monterey, California. – p.6/34
 Stiff ODEs
Explicit methods for solving ordinary differential equations had particular
difficulty when the problem was stiff.
The use of implicit methods based on backward difference formula was
brought to the notice of linear algebraists by the landmark talk by Bill Gear
at the IFIP Conference in Edinburgh in 1968.
The equations that were solved in Gear’s backward difference formula
were of the form

(I − αhJ)x = b
where J was the Jacobian, h the stepsize, and α a parameter.
These equations could be very evil but very sparse.
Liniger at IBM Yorktown and Curtis at Harwell both had equations of this
kind to solve.

The SIAM Conference on Applied Linear Algebra. October 26-29, 2009. Monterey, California. – p.7/34
 Linear Programming

The linear programming community was one of the first to embrace

sparsity
1947 The simplex method. George Dantzig
1954 The product form of the inverse. Dantzig and Orchard-Hays
1955-56 Orchard-Hays implements sparse LP codes on the Johnniac at
RAND
1957 The elimination form of the inverse. Harry Markowitz
1961 Martin Beale doing sparse programming at CEIR (UK) later
Scientific Control (SCICON)

The SIAM Conference on Applied Linear Algebra. October 26-29, 2009. Monterey, California. – p.8/34
Ancient times

Harry M. Markowitz
The SIAM Conference on Applied Linear Algebra. October 26-29, 2009. Monterey, California. – p.9/34
The Birth

The SIAM Conference on Applied Linear Algebra. October 26-29, 2009. Monterey, California. – p.10/34
The Birth

The SIAM Conference on Applied Linear Algebra. October 26-29, 2009. Monterey, California. – p.10/34
 The Birth
The first Sparse Matrix Symposium was organized at IBM Yorktown
Heights in 1968 by staff of the Mathematical Sciences Department: Robert
A. Brayton, Fred G. Gustavson, Alan J. Hoffman, and Philip S. Wolfe, and
Ralph A. Willoughby, who edited the Proceedings.

I have always thought of Ralph as the “father” of sparse matrix technology.

Obituary written by Jane Cullum in SIAM News 35 (7), September 2003.

The SIAM Conference on Applied Linear Algebra. October 26-29, 2009. Monterey, California. – p.11/34
 The Birth

Ralph Willoughby 1923 - 2001 [courtesy of family and SIAM]

The SIAM Conference on Applied Linear Algebra. October 26-29, 2009. Monterey, California. – p.12/34
The Birth

The SIAM Conference on Applied Linear Algebra. October 26-29, 2009. Monterey, California. – p.13/34
The Birth

The SIAM Conference on Applied Linear Algebra. October 26-29, 2009. Monterey, California. – p.14/34
The Birth

The SIAM Conference on Applied Linear Algebra. October 26-29, 2009. Monterey, California. – p.15/34
The Birth

The SIAM Conference on Applied Linear Algebra. October 26-29, 2009. Monterey, California. – p.16/34
The Birth

The SIAM Conference on Applied Linear Algebra. October 26-29, 2009. Monterey, California. – p.17/34
 The Birth
A number of “currently hot" topics were discussed at the 1968 meeting!
Hybrid methods
MC64 scaling before Olschowka and Neumaier (1996)
Use of single precision arithmetic to get full accuracy
Permutations to block triangular and bordered block diagonal form and
“tearing” in general
Modified Markowitz

The SIAM Conference on Applied Linear Algebra. October 26-29, 2009. Monterey, California. – p.18/34
 The Birth
At this time perhaps the only “general purpose” code was that of Fred
Gustavson at IBM Yorktown who used a generated code approach for
solving problems coming from stiff ODEs. This code was later called
GNSOIN.
Papers by Brayton, Gustavson, and Willoughby (1970), Gustavson,
Liniger, and Willoughby on GNSOIN (1970), and by Hachtel, Brayton,
and Gustavson on variability typing (1971).
Shortly afterwards, Curtis and Reid developed codes MA17 (symmetric
positive definite) and MA18 (unsymmetric), the latter for the same
application as above.

The SIAM Conference on Applied Linear Algebra. October 26-29, 2009. Monterey, California. – p.19/34
The Christening

The SIAM Conference on Applied Linear Algebra. October 26-29, 2009. Monterey, California. – p.20/34
The Christening

The SIAM Conference on Applied Linear Algebra. October 26-29, 2009. Monterey, California. – p.20/34
 The infancy

Around this time, the first theses on sparse direct methods appeared.

Donald Rose 1970 Harvard Symmetric elimination on sparse

positive definite systems and the po-
tential flow network problem
Alan George 1971 Stanford Computer implementation of the
finite-element method
Iain Duff 1972 Oxford Analysis of sparse systems
Andrew Sherman 1975 Yale On the efficient solution of sparse
systems of linear and non-linear
equations

The SIAM Conference on Applied Linear Algebra. October 26-29, 2009. Monterey, California. – p.21/34
 The infancy

By the mid 1970s, there were only really three widely available sparse
matrix packages: HSL, SPARSPAK, and YSMP.
Harwell Subroutine Library (HSL)
Main reason for developing these codes was to enable Alan Curtis to solve
stiff ODE problems. First codes were MA17 (symmetric positive definite)
and MA18 (unsymmetric) by Curtis and Reid in 1971 and then MA28
(unsymmetric) by Duff in 1977 and MA27 (symmetric) by Duff and Reid
in 1982.
YSMP (1975) LLL User’s Guide by Andrew Sherman. Paper by Eisenstat,
Gursky, Schultz, and Sherman in 1982
SPARSPAK (1978) University of Waterloo. George and Liu (then Ng).
By 1982 a sparse software catalogue also lists codes from IBM, Bell Labs,
CRAY, FPS, SSLEST and Y12M from DTU, COSMIC, NSPFAC, ...

The SIAM Conference on Applied Linear Algebra. October 26-29, 2009. Monterey, California. – p.22/34
 Maturity

SPARSE BOOKS
1973 Tewarson Sparse Matrices
1976 Brameller, Allan and Hamam Sparsity
1981 George and Liu Computer Solution of Large Sparse
Positive Definite Systems
1983 Østerby and Zlatev Direct Methods for Sparse Matrices
1984 Pissanetsky Sparse Matrix Technology
1986 Duff, Erisman and Reid Direct Methods for Sparse Matrices
1991 Zlatev Computational Methods for General
Sparse Matrices
2006 Davis Direct Methods for Sparse Linear
Systems

The SIAM Conference on Applied Linear Algebra. October 26-29, 2009. Monterey, California. – p.23/34
 Maturity

Conference Proceedings
1969 Willoughby Sparse Matrix Proceedings
1971 Reid Large Sparse Sets of Linear Equations
1972 Rose and Willoughby Sparse Matrices and their Applica-
tions
1973 Himmelblau Decomposition of Large-Scale Prob-
lems
1976 Bunch and Rose Sparse Matrix Computations
1977 Barker Sparse Matrix Techniques
1979 Duff and Stewart Sparse Matrix Proceedings 1978
1981 Duff Sparse Matrices and their Uses
1985 Evans Sparsity and its Applications

The SIAM Conference on Applied Linear Algebra. October 26-29, 2009. Monterey, California. – p.24/34
 Maturity

Other publications and meetings that established the field.

1975 Symposium on Sparse Matrix Computations. Argonne
1977 A Survey of Sparse Matrix Research. Proc IEEE. Duff
1980 Direct Methods for solving large sparse systems. SIAM News.
George
1982 Sparse Matrix Symposium 1982. Fairfield Glade, Tennessee.
1989 SIAM Symposium on Sparse Matrices. Salishan, Oregon
1990 Sparsity in Large Scientific Computations. IBM Europe Institute.
Oberlech, Austria.
1995 International Linear Algebra Year at CERFACS. Direct methods
workshop 1995. BIT 37, 1997.
1996 SIAM Meeting on Sparse Matrices. Coeur d’Alene, Idaho.
A number of CSC (Combinatorial Scientific Computing) Meetings.
Roughly biennial.
The SIAM Conference on Applied Linear Algebra. October 26-29, 2009. Monterey, California. – p.25/34
 Turbulent times
Gordon Research Conference on Flow in Permeable Media. Proctor
Academy, Andover, NH. July 30 - Aug 3, 1984

The SIAM Conference on Applied Linear Algebra. October 26-29, 2009. Monterey, California. – p.26/34
 Turbulent times
Gordon Research Conference on Flow in Permeable Media. Proctor
Academy, Andover, NH. July 30 - Aug 3, 1984

Talk of ISD
The solution of large sparse asymmetric matrices by direct methods.

The SIAM Conference on Applied Linear Algebra. October 26-29, 2009. Monterey, California. – p.26/34
 Turbulent times
Gordon Research Conference on Flow in Permeable Media. Proctor
Academy, Andover, NH. July 30 - Aug 3, 1984

Talk of ISD
The solution of large sparse asymmetric matrices by direct methods.

Punchline
Real men use direct methods

The SIAM Conference on Applied Linear Algebra. October 26-29, 2009. Monterey, California. – p.26/34
 Maturity

Grid dimensions Matrix order Work to factorize Factor storage

k×k k2 k3 k 2 log k

k×k×k k3 k6 k4

O complexity of direct method on 2D and 3D grids.

The SIAM Conference on Applied Linear Algebra. October 26-29, 2009. Monterey, California. – p.27/34
Maturity ... multifrontal

The SIAM Conference on Applied Linear Algebra. October 26-29, 2009. Monterey, California. – p.28/34
 Speelpenning’s contribution .. 1973

Background is very much in frontal methods (Irons 1970) and

finite-element computations
Makes break from wavefront/frontal ordering
Uses generalized element formulation
Combines elements on elimination to form larger elements
Does multiple eliminations within a front
No numerical experiments
Preamble to 1978 reprint of 1973 report references use in DIANA
system at TNO-IBBC for solving systems from structural analysis of
order 5000 with “bandwidth” about 400

The SIAM Conference on Applied Linear Algebra. October 26-29, 2009. Monterey, California. – p.29/34
 Early multifrontal developments

Sherman (1975) : compressed storage scheme

Eisenstat, Schultz, and Sherman (1976) : element merge tree
Peters (1980) : substructuring (claims not generalized element)
George and Liu (1980) : generalized element approach for minimum
degree
Duff and Reid (1983) : paper on background to multifrontal code MA27
Reid (1984) : TREESOLV (out-of-core)
Duff (1986) : prototype parallel implementation on Alliant FX/8
Liu (1992) : detailed description of multifrontal method

The SIAM Conference on Applied Linear Algebra. October 26-29, 2009. Monterey, California. – p.30/34
 MA27 ... 1982
Main contributions
Full implementation of tree for all phases (including node
amalgamation)
Use of stack for intermediate frontal matrices
Symmetric indefinite factorization using Duff, Reid, Munksgaard, and
Nielsen’s (1979) sparse adaptation of a variant of two-by-two pivoting
strategy of Bunch and Parlett (1971). This was perhaps the first
implementation of rook pivoting that Saunders used later in LUSOL.
Exploitation of vector machines (pre-dated higher-level BLAS) using
experience from frontal solvers (MA32)

Note that the only codes with which we could then compare MA27 were
MA17, YSMP and SPARSPAK (plus a code of Munksgaard)

The SIAM Conference on Applied Linear Algebra. October 26-29, 2009. Monterey, California. – p.31/34
Senility?

The SIAM Conference on Applied Linear Algebra. October 26-29, 2009. Monterey, California. – p.32/34
 Senility?

HYBRID METHODS
COMBINING DIRECT AND ITERATIVE METHODS
(can be thought of as sophisticated preconditioning)
Multigrid
Using direct method as coarse grid solver.
Domain Decomposition
Using direct method on local subdomains and “direct” preconditioner on
interface.
Block Iterative Methods

Direct solver on sub-blocks.

Partial factorization as preconditioner

Factorization of nearby problem as a preconditioner

The SIAM Conference on Applied Linear Algebra. October 26-29, 2009. Monterey, California. – p.32/34
 Conclusions

Sparse Direct Methods have had a short, distinguished,

and very active history
There are still many future research challenges

A range of techniques involving both sparse direct and a

range of sparse iterative solvers is required to solve the
really challenging problems of the future

The SIAM Conference on Applied Linear Algebra. October 26-29, 2009. Monterey, California. – p.33/34
Conclusions

THANK YOU
for your attention

The SIAM Conference on Applied Linear Algebra. October 26-29, 2009. Monterey, California. – p.34/34