Symmetry reduces the amount of work required to decompose symmetric systems into triangular factors. It does not reduce the work required to actually solve the system from an existing triangular factorization. Implementing forward substitution for a symmetric decomposition boils down to making sure implicit data (i.e. the portion of the symmetric factorization that is not not physically stored) is correctly derived from the explicitly stored data. See Section 7.2 for a discussion of implicit data in symmetric LU decompositions.
Forward substitution solves lower triangular systems. When L is available, the symmetric and asymmetric solution algorithms are identical. See Section 5.1 for the general implementation of forward substitution. If L is stored in a linear array, use Equation 86 or Equation 87 for indexing.
The case in which U is the explicit factorization is examined more closely. If U is an n × n matrix containing the upper triangular factors of asymmetric matrix, then L is unit lower triangular and obtained from U via Equation 81 with lij being one. Substituting Equation 81 into Equation 62 yields:
| $$y_{i}=b_{i}-\displaystyle\sum_{j=i+1}^{n}\frac{u_{ji}y_{j}}{u_{jj}},\mbox{ where }1 \leq i \leq n$$ | (88) |
Algorithm 14 implements Equation 88, i.e. the inner product formulation of forward substitution for symmetric systems whose upper triangular factors are available.
| for $i=1,\cdots,n$ |
| $\alpha=b_{i}$ |
| for $j=1,\cdots,i-1$ |
| $\alpha=\alpha-\displaystyle\frac{u_{ji}}{u_{jj}}y_{j}$ |
| $y_{i}=\alpha$ |
If U is stored in alinear array with zero based indexing, the inner product formulation of forward substitution is implemented by Algorithm 15.
| for $i=0,\cdots,n-1$ |
| $\alpha=\mathbf{b}{[i]}$ |
| for $j=0,\cdots,i-1$ |
| $\alpha=\alpha-\displaystyle\frac{\mathbf{u}{[jn-j(j+1)/2+i]}}{\mathbf{u}% {[jn-j(j+1)/2+j]}}\cdot\mathbf{y}{[j]}$ |
| $\mathbf{y}{[i]}=\alpha$ |
Algorithm 16 implements the outer product formulation of forward substitution for symmetric systems whose upper triangular factors are available.
| for $i=1,\cdots,n$ |
| $y_{i}=b_{i}$ |
| $\alpha=\displaystyle\frac{y_{i}}{u_{ii}}$ |
| for $k=i+1,\cdots,n$ |
| $b_{k}=b_{k}-\alpha u_{ik}$ |
This differs from the asymmetric implementation in that
| $$l_{ii}=1\operatorname{and}l_{ki}=\frac{u_{ik}}{u_{ii}},\mbox{ where }i \neq k$$ | (89) |
Therefore, the initial division by lii is omitted and the division by uii is pulled out of the k loop. The outer product formulation of forward substitution where U is stored in a linear array with zero based indexing is realized by Algorithm 17.
| for $i=0,\cdots,n-1$ |
| $\mathbf{y}[i] = \mathbf{b}[i]$ |
| $\alpha=\displaystyle\frac{\mathbf{y}{[i]}}{\mathbf{u}{[in-i(i+1)/2+i% ]}}$ |
| for $k=i+1,\cdots,n-1$ |
| $\mathbf{b}[k] = \mathbf{b}[k]-\alpha\cdot\mathbf{u}[in-i(i+1)/2+k]$ |