7.5 Crout’s Method for Symmetric Matrices

If A is a symmetric n × n matrix, Algorithm 12 computes – one column at a time – the lower triangular matrix L that results from a Crout decomposition. The lower triangular portion of A is overwritten by L.

Algorithm 12: Crout’s Method - Symmetric Implementation
for j=1,,nj=1,\cdots,n
   for i=1,,j-1i=1,\cdots,j-1
      wi=ajiaiiw_{i}=\displaystyle\frac{a_{ji}}{a_{ii}}
   for i=j,,ni=j,\cdots,n
      α=aij\alpha=a_{ij}
      for p=1,,j-1p=1,\cdots,j-1
         α=α-aipwp\alpha=\alpha-a_{ip}w_{p}
      aij=αa_{ij}=\alpha

The algorithm uses a working vector w of length n to construct the relevant portion of column j from U at each stage of the factorization. Elements from U are derived using Equation 63.

If the lower triangular portion of A is stored in a linear array, Algorithm 13 results in the same factorization (assuming zero based subscripting). The algorithm overwrites A with L.

Algorithm 13: Crout’s Method - Symmetric, Array Based
for j=0,,n-1j=0,\cdots,n-1
   for i=0,,j-1i=0,\cdots,j-1
      𝐰[j]=𝐚[jn-j(j+1)/2+i]𝐚[jn-j(j+1)/2+j]\mathbf{w}\text{[j]}=\displaystyle\frac{\mathbf{a}\text{[jn-j(j+1)/2+i]}}{% \mathbf{a}\text{[jn-j(j+1)/2+j]}}
   for i=j,,n-1i=j,\cdots,n-1
      α=𝐚[i(i+1)/2+j]\alpha=\mathbf{a}\text{[i(i+1)/2+j]}
      for p=1,,i-1p=1,\cdots,i-1
         α=α\alpha=\alpha - 𝐚\mathbf{a}[i(i+1)/2+p] \cdot 𝐰\mathbf{w}[p]
      𝐚\mathbf{a}[i(i+1)/2+j] = α\alpha

For the general implementation of Crout’s method, see Section 4.3.