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 70 or Equation 71 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 a symmetric matrix, then L is unit lower triangular and obtained from U via Equation 65 with ljj{}_{jj} being one. Substituting Equation 65 into Equation 46 yields:

yi=bi-j=i+1nujiyjujj,where1iny_{i}=b_{i}-\displaystyle\sum_{j=i+1}^{n}\frac{u_{ji}y_{j}}{u_{jj}},% \operatorname{where}1\leq i\leq n (72)

Algorithm 14 implements Equation 72, 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,,ni=1,\cdots,n
   α=bi\alpha=b_{i}
   for j=1,,i-1j=1,\cdots,i-1
      α=α-ujiujjyj\alpha=\alpha-\displaystyle\frac{u_{ji}}{u_{jj}}y_{j}
   yi=αy_{i}=\alpha

If U is stored in a linear 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,,n-1i=0,\cdots,n-1
   α=𝐛[i]\alpha=\mathbf{b}\text{[i]}
   for j=0,,i-1j=0,\cdots,i-1
      α=α-𝐮[jn-j(j+1)/2+i]𝐮[jn-j(j+1)/2+j]𝐲[j]\alpha=\alpha-\displaystyle\frac{\mathbf{u}\text{[jn-j(j+1)/2+i]}}{\mathbf{u}% \text{[jn-j(j+1)/2+j]}}\cdot\mathbf{y}\text{[j]}
   𝐲[i]=α\mathbf{y}\text{[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,,ni=1,\cdots,n
   yi=biy_{i}=b_{i}
   α=yiuii\alpha=\displaystyle\frac{y_{i}}{u_{ii}}
   for k=i+1,,nk=i+1,\cdots,n
      bk=bk-αuikb_{k}=b_{k}-\alpha u_{ik}

This differs from the asymmetric implementation in that

lii=1andlki=uikuii,whereikl_{ii}=1\operatorname{and}l_{ki}=\frac{u_{ik}}{u_{ii}},\operatorname{where}i\neq k (73)

Therefore, the initial division by lii{}_{ii} is omitted and the division by uii{}_{ii} 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,,n-1i=0,\cdots,n-1
   𝐲\mathbf{y}[i] = 𝐛\mathbf{b}[i]
   α=𝐲[i]𝐮[in-i(i+1)/2+i]\alpha=\displaystyle\frac{\mathbf{y}\text{[i]}}{\mathbf{u}\text{[in-i(i+1)/2+i% ]}}
   for k=i+1,,n-1k=i+1,\cdots,n-1
      𝐛\mathbf{b}[k] = 𝐛\mathbf{b}[k] -α𝐮-\alpha\cdot\mathbf{u}[in-i(i+1)/2+k]