51
Recursive approach in sparse matrix LU
factorization
Jack Dongarra, Victor Eijkhout and
Piotr Łuszczek
∗
University of Tennessee, Department of Computer
Science, Knoxville, TN 37996-3450, USA
Tel.: +865 974 8295; Fax: +865 974 8296
This paper describes a recursive method for the LU factoriza-
tion of sparse matrices. The recursive formulation of com-
mon linear algebra codes has been proven very successful in
dense matrix computations. An extension of the recursive
technique for sparse matrices is presented. Performance re-
sults given here show that the recursive approach may per-
form comparable to leading software packages for sparse ma-
trix factorization in terms of execution time, memory usage,
and error estimates of the solution.
1. Introduction
Typically, a system of linear equations has the form:
Ax = b, (1)
where A is n by n real matrix (A ∈ R
n×n
), and x
and b are n-dimensional real vectors (b, x ∈ R
n
). The
values of A and b are known and the task is to find
x satisfying Eq. (1). In this paper, it is assumed that
the matrix A is large (of order commonly exceeding
ten thousand) and sparse (there are enough zero entries
in A that it is beneficial to use special computational
methods to factor the matrix rather than to use a dense
code). There are two common approaches that are used
to deal with such a case, namely, iterative [33] and
direct methods [17].
Iterative methods, in particular Krylov subspace
techniques such as the Conjugate Gradient algorithm,
are the methods of choice for the discretizations of el-
liptic or parabolic partial differential equations where
∗
Corresponding author: Piotr Luszczek, Department of Computer
Science, 1122 Volunteer Blvd., Suite 203, Knoxville, TN 37996-
3450, USA. Tel.: +865 974 8295; Fax: +865 974 8296; E-mail:
luszczek@cs.utk.edu.
the resulting matrix is often guaranteed to be positive
definite or close to it. However, when the linear sys-
tem matrix is strongly unsymmetric or indefinite, as
is the case with matrices originating from systems of
ordinary differential equations or the indefinite matri-
ces arising from shift-invert techniques in eigenvalue
methods, one has to revert to direct methods which are
the focus of this paper.
In direct methods, Gaussian elimination with partial
pivoting is performedto find a solution of Eq. (1). Most
commonly, the factored form of A is given by means
of matrices L, U, P and Q such that:
LU = PAQ, (2)
where:
– L is a lower triangular matrix with unitary diago-
nal,
– U is an upper triangular matrix with arbitrary di-
agonal,
– P and Q are row and column permutation matri-
ces, respectively (each row and column of these
matrices contains single a non-zero entry which is
1, and the following holds: PP
T
= QQ
T
= I,
with I being the identity matrix).
A simple transformation of Eq. (1) yields:
(PAQ)Q
−1
x = Pb, (3)
which in turn, after applying Eq. (2), gives:
LU(Q
−1
x)=Pb, (4)
SolutiontoEq.(1)maynowbeobtainedintwosteps:
Ly = Pb (5)
U(Q
−1
x)=y (6)
and these steps are performed through forward/back-
ward substitution since the matrices involved are trian-
gular. The most computationally intensive part of solv-
ing Eq. (1) is the LU factorization defined by Eq. (2).
This operation has computational complexity of or-
der O(n
3
) when A is a dense matrix, as compared
to O(n
2
) for the solution phase. Therefore, optimiza-
Scientific Programming 9 (2001) 51–60
ISSN 1058-9244 / $8.00
2001, IOS Press. All rights reserved
评论0