8 Sparse Matrices 8.5 Creating PAP^T from a Symbolic Factorization 8.7 Solving Sparse Linear Systems

8.6 Numeric Factorization of Sparse Matrices

Numeric factorization algorithms work with nonzero values of a sparse matrix A and the data structures resulting from symbolic factorization to compute the factorization

\mathbf{LU=PAP^{T}}

Algorithms discussed in the current section act on a sparse n × n matrix A. They assume that

A already reflects the pivot strategy defined by P and P^T, i.e. the algorithms pivot down the diagonal.
A has zero-valued entries at fill-up locations.
A is maintained in a sparse data structure supporting the get, scan, and put operators. Communication with the data structure is maintained through the buffer a in the normal manner.

Techniques for creating the requisite pivot ordered, fill-up augmented A matrix (i.e. PAP^T) are discussed in Section 8.5.

Algorithm 25 uses Doolittle’s method to compute the LU decomposition of a sparse matrix A.The algorithm overwrites A with LU.

Algorithm 25: LU Decomposition of a Sparse Matrix by Doolittle’s Method ^†

for

i=2,\cdots,n

while[

j=

scan(

\mathbf{A}

,row

i,1,n,a

)]

\neq

finished

push_scan

\alpha=a

m=\text{ }

min(

i-1,j-1

)

while[

p=

scan(

\mathbf{A}

,row

i,1,m,a

)]

\neq

finished

b=a

if get(

\mathbf{A},p,j,a

)is successful

\alpha=\alpha-ab

i>j

get(

\mathbf{A},j,j,a

)

\alpha=\displaystyle\frac{\alpha}{a}

a=\alpha

put(

\mathbf{A},i,j,a

)

pop_scan

^†Operations and data types are defined in Section 8.2.

Algorithm 26 uses Doolittle’s method to compute U, the upper triangular factors, of a symmetric sparse matrix A. It is assumed that A is initially stored as an upper triangular matrix. The algorithm overwrites A with U. The vector w is used to construct the nonzero entries of each column from U. The vector c contains cursors to the row in L with which the entries of A are associated, e.g. if w_k contains l_ji then c_k is j.

Algorithm 26: LU Decomposition of Sparse Symmetric Matrix by Doolittle’s Method ^†

for

i=1,\cdots,n

k=0

for

j=1,\cdots,i-1

if get(

\mathbf{A},j,i,a

)is successful

k=k+1

w_{k}=a

c_{k}=j

get(

\mathbf{A},j,j,a

)

w_{k}=\displaystyle\frac{w_{k}}{a}

while[

j=

scan(

\mathbf{A}

,row

i,i,n,a

)]

\neq

finished

\alpha=a

for

p=1,\cdots,k

if get(

\mathbf{A},c_{p},j,a

)is successful

\alpha=\alpha-aw_{p}

a=\alpha

put(

\mathbf{A},i,j,a

)

^†Operations and data types are defined in Section 8.2.

Doolittle’s method for full, asymmetric matrices, is examined in Section 4.2. Doolittle’s method for full, symmetric matrices is discussed in Section 7.4.