7.4 Doolittle’s Method for Symmetric Matrices

If A is a symmetric n × n matrix, Algorithm 10 computes – one row at a time – the upper triangular matrix U that results from a Doolittle decomposition. The upper triangular portion of A is overwritten by U.

Algorithm 10: Doolittle’s Method - Symmetric Implementation
for i=1,,ni=1,\cdots,n
   for j=1,,i-1j=1,\cdots,i-1
      wj=ajiajjw_{j}=\displaystyle\frac{a_{ji}}{a_{jj}}
   for j=i,,nj=i,\cdots,n
      α=aij\alpha=a_{ij}
      for p=1,,i-1p=1,\cdots,i-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 row i from L at each stage of the factorization. Elements from L are derived using Equation 65.

If the upper triangular portion of A is stored in a linear array, Algorithm 11 results in the same factorization (assuming zero based subscripting). The upper triangular portion ofA is overwritten by U.

Algorithm 11: Doolittle’s Method - Symmetric, Array Based
for i=0,,n-1i=0,\cdots,n-1
   for j=0,,i-1j=0,\cdots,i-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 j=i,,n-1j=i,\cdots,n-1
      α=𝐚\alpha=\mathbf{a}[in-i(i+1)/2+j]
      for p=1,,i-1p=1,\cdots,i-1
         α=α-𝐚\alpha=\alpha-\mathbf{a}[in-i(i+1)/2+p]𝐰\cdot\mathbf{w}[p]
      𝐚\mathbf{a}[in-i(i+1)/2+j] = α\alpha

For the general implementation of Doolittle’s method, see Section 4.2.