8.7 Solving Sparse Linear Systems

As we have seen, the LU decomposition of the matrix PAQ is used to solve the linear system of equations Ax = b by sequentially solving Equation 42 through Equation 45, which are repeated below.

𝐲\displaystyle\mathbf{y} =𝐏𝐛\displaystyle=\mathbf{Pb}
𝐋𝐜\displaystyle\mathbf{Lc} =𝐲\displaystyle=\mathbf{y}
𝐔𝐳\displaystyle\mathbf{Uz} =𝐜\displaystyle=\mathbf{c}
𝐱\displaystyle\mathbf{x} =𝐐𝐳\displaystyle=\mathbf{Qz}

Sparsity techniques benefit this process. The algorithms presented in this section assume:

  • An LU decomposition of PAQ exists.

  • The factorization is maintained in a sparse matrix data structure which supports the scan and get operators. These operators are described in Section 8.2.1.

  • The row permutations P are represented by a mapping ψ\psi whose domain is a row index in A and whose range is the corresponding row index in PAQ.

  • The column permutations Q are represented by a mapping τ\tau whose domain is a column index in PAQ and whose range is the corresponding column index in A.

For example, Section 8.6 describes algorithms that create numeric factorizations satisfying the first two of these assumptions. Section 8.4.1 describes an algorithm for obtaining the row and column permutations corresponding to a minimum degree pivot strategy.

For simplicity of exposition, it is assumed that the vectors b, c, x, y, and z are stored in linear arrays. However, any data structure that lets you retrieve and update an element of a vector based on its index will suffice.

8.7.1 Permute the Constant Vector

The equation y = Pb is efficiently implemented using the mapping ψ\psi to permute the elements of b. Algorithm 27 illustrates this procedure.

Algorithm 27: Permute b to order P
for i=1,,ni=1,\cdots,n
   k=ψ(i)k=\psi(i)
   yi=bky_{i}=b_{k}

8.7.2 Sparse Forward Substitution

The lower triangular system of equations Lc = y is solved by forward substitution. Algorithm 28 implements the inner product formulation of forward substitution on a sparse L and full vectors c and y.

Algorithm 28: Sparse Forward Substitution
for i=1,,ni=1,\cdots,n
   α=yi\alpha=y_{i}
   while[j=j=scan(𝐋\mathbf{L},row i,1,i-1,li,1,i-1,l)] \neq finished
      α\alpha =  α-lyi\alpha-ly_{i}
   get(𝐋,i,i,l\mathbf{L},i,i,l)
   ci=αlc_{i}=\displaystyle\frac{\alpha}{l}

The sequencing of the operations permits the use of a single vector y. Operations that update ci would overwrite yi instead. If L is unit lower triangular, division by the diagonal element lii is omitted.

Algorithm 29 implements an outer product formulation of forward substitution for use with symmetric systems whose upper triangular factors are available. See Section 5.3 and Section 7.6 for additional information.

Algorithm 29: Symmetric Sparse Forward Substitution
for i=n,,1i=n,\cdots,1
   ci=yic_{i}=y_{i}
   get(𝐔,i,i,u\mathbf{U},i,i,u)
   α=ciu\alpha=\displaystyle\frac{c_{i}}{u}
   while[k=k=scan(𝐔\mathbf{U},row i,i+1,n,ui,i+1,n,u)] \neq finished
      yk=yk-uαy_{k}=y_{k}-u\alpha

8.7.3 Sparse Backward Substitution

The upper triangular system of equations Uz = c is solved by backward substitution. Algorithm 30 implements the inner product formulation of backward substitution on a sparse U and full vectors c and z.

Algorithm 30: Sparse Back Substitution
for i=n,,1i=n,\cdots,1
   α=ci\alpha=c_{i}
   while[j=j=scan(𝐔\mathbf{U},row i,i+1,n,ui,i+1,n,u)] \neq finished
      α=α-uzj\alpha=\alpha-uz_{j}
   get(𝐔,i,i,u\mathbf{U},i,i,u)
   zi=αuz_{i}=\displaystyle\frac{\alpha}{u}

The sequencing of the operations permits the use of a single vector c. Operations that update zi would overwrite ci instead. If U is unit upper triangular, division by the diagonal element uii is omitted.

8.7.4 Permute the Solution Vector

The equation x = Qz is efficiently implemented using the mapping τ\tau to permute the elements of z. Algorithm 31 illustrates this process.

Algorithm 31: Permute x to order Q
for i=1,,ni=1,\cdots,n
   k=τ(i)k=\tau(i)
   xi=zkx_{i}=z_{k}