7.7 Backward Substitution for Symmetric Systems

As stated previously, 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 backward substitution for a symmetric decomposition reduces to making sure that 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.7.1 Back Substitution Using Upper Triangular Factors

Backward substitution solves upper triangular systems. When U is available, the symmetric and asymmetric solution algorithms are identical. See Section 5.2 for the general implementation of backward substitution. If U is stored in a linear array, use Equation 70 or Equation 71 for indexing.

7.7.2 Back Substitution Using Lower Triangular Factors

The case in which L is the explicit factorization merits further attention. If L is an n × n matrix containing the upper triangular factors of a symmetric matrix, then U is unit lower triangular and obtained from L via Equation 63 with uii{}_{ii} being one. Substituting Equation 63 into Equation 47 yields:

xi=yi-j=i+1nljixjljj,wherei=n,n-1,,1x_{i}=y_{i}-\displaystyle\sum_{j=i+1}^{n}\frac{l_{ji}x_{j}}{l_{jj}},% \operatorname{where}i=n,n-1,\cdots,1 (74)

Algorithm 18 implements this inner product formulation of backward substitution for symmetric systems whose lower triangular factors are available.

Algorithm 18: Symmetric Back Substitution using Lower Triangular Factors
for i=1,,ni=1,\cdots,n
   for j=1,,i-1j=1,\cdots,i-1

If L is stored in a linear array with zero based indexing, the inner product formulation of back substitution is stated in Algorithm 19.

Algorithm 19: Symmetric Backward Substitution using L with Array Storage
for i=n-1,,0i=n-1,\cdots,0
   for j=i+1,,n-1j=i+1,\cdots,n-1
      α=α-𝐥[jn-j(j+1)/2+i]𝐥[jn-j(j+1)/2+j])𝐱[j]\alpha=\alpha-\displaystyle\frac{\mathbf{l}\text{[jn-j(j+1)/2+i]}}{\mathbf{l}% \text{[jn-j(j+1)/2+j])}}\cdot\mathbf{x}\text{[j]}