As we have seen, the LU decomposition of the matrix PAQ isused to solve the linear system of equations Ax = b by sequentially solving Equation 58 through Equation 61, which are repeated below.
| $$\mathbf{y}=\mathbf{Pb}$$ $$\mathbf{Lc}=\mathbf{y}$$ $$\mathbf{Uz}=\mathbf{c}$$ $$\mathbf{x}=\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.
The equation y = Pb is efficiently implemented using the mapping $\psi$ to permute the elements of b. Algorithm 27 illustrates this procedure.
| for $i=1,\cdots,n$ |
| $k=\psi(i)$ |
| $y_{i}=b_{k}$ |
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.
| for $i=1,\cdots,n$ |
| $\alpha=y_{i}$ |
| while [ $j=$ scan( $\mathbf{L}$,row $i,1,i-1,l$ ) ] $\neq$ finished |
| $\alpha=\alpha-ly_{i}$ |
| get( $\mathbf{L},i,i,l$ ) |
| $c_{i}=\displaystyle\frac{\alpha}{l}$ |
| $^\dagger$ Operations and data types are defined in Section 8.2.1. |
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.
| for $i=n,\cdots,1$ |
| $c_{i}=y_{i}$ |
| get( $\mathbf{U},i,i,u$ ) |
| $\alpha=\displaystyle\frac{c_{i}}{u}$ |
| while [ $k=$ scan( $\mathbf{U}$,row $i,i+1,n,u$ ) ] $\neq$ finished |
| $y_{k}=y_{k}-u\alpha$ |
| $^\dagger$ Operations and data types are defined in Section 8.2.1. |
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.
| for $i=n,\cdots,1$ |
| $ \alpha=c_{i}$ |
| while [ $j=$scan( $\mathbf{U}$,row $i,i+1,n,u$ ) ] $\neq$ finished |
| $\alpha=\alpha-uz_{j}$ |
| get( $\mathbf{U},i,i,u$ ) |
| $z_{i}=\displaystyle\frac{\alpha}{u}$ |
| $^\dagger$ Operations and data types are defined in Section 8.2.1. |
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.
The equation x = Qz is efficiently implemented using the mapping $\tau$ topermute the elements of z. Algorithm 31 illustrates this process.
| for $i=1,\cdots,n$ |
| $k=\tau(i)$ |
| $x_{i}=z_{k}$ |