7.6 Forward Substitution for Symmetric Systems

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.

7.6.1 Forward Substitution Using Lower Triangular Factors

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.

7.6.2 Forward Substitution Using Upper Triangular Factors

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.

Algorithm 14: Symmetric Forward Substitution via Upper Triangular Factors
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.

Algorithm 15: Symmetric Forward Substitution using U with Array Storage
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.

Algorithm 16: Symmetric Forward Substitution using U, Outer Product
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.

Algorithm 17: Symmetric Forward Substitution using U, Outer Product, Array
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]$