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 PT, 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. PAPT) 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,,ni=2,\cdots,n
   while[j=j=scan(𝐀\mathbf{A},row i,1,n,ai,1,n,a)] \neq finished
      push_scan
      α=a\alpha=a
      m= m=\text{ }min(i-1,j-1i-1,j-1)
      while[p=p=scan(𝐀\mathbf{A},row i,1,m,ai,1,m,a)] \neq finished
          b=ab=a
          if get(𝐀,p,j,a\mathbf{A},p,j,a)is successful
              α=α-ab\alpha=\alpha-ab
      if i>ji>j
          get(𝐀,j,j,a\mathbf{A},j,j,a)
          α=αa\alpha=\displaystyle\frac{\alpha}{a}
      a=αa=\alpha
      put(𝐀,i,j,a\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 wk contains lji then ck is j.

Algorithm 26: LU Decomposition of Sparse Symmetric Matrix by Doolittle’s Method 
for i=1,,ni=1,\cdots,n
   k=0k=0
   for j=1,,i-1j=1,\cdots,i-1
      if get(𝐀,j,i,a\mathbf{A},j,i,a)is successful
          k=k+1k=k+1
          wk=aw_{k}=a
          ck=jc_{k}=j
          get(𝐀,j,j,a\mathbf{A},j,j,a)
          wk=wkaw_{k}=\displaystyle\frac{w_{k}}{a}
   while[j=j=scan(𝐀\mathbf{A},row i,i,n,ai,i,n,a)] \neq finished
       α=a\alpha=a
       for p=1,,kp=1,\cdots,k
          if get(𝐀,cp,j,a\mathbf{A},c_{p},j,a)is successful
             α=α-awp\alpha=\alpha-aw_{p}
       a=αa=\alpha
       put(𝐀,i,j,a\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.