An LU decomposition of A may be obtained by applying the definition of matrix multiplication to the equation A = LU. If L is unit lower triangular and U is upper triangular, then
Rearranging the terms of Equation 29 yields
Jointly Equation 30 and Equation 31 are referred to as Doolittle’s method of computing the LU decomposition of A. Algorithm 2 implements Doolittle’s method. Calculations are sequenced to compute one row of L followed by the corresponding row of U until A is exhausted. You will observe that the procedure overwrites the subdiagonal portions of A with L and the upper triangular portions of A with U.
Figure 1 depicts the computational sequence associated with Doolittle’s method.
Two loops in the Doolittle algorithm are of the form
These loops determine an element of the factorization lij or uij by computing the dot product of a partial column vector in U and partial row vector in L. As such, the loops perform an inner product accumulation. These computations have a numerical advantage over the gradual accumulation of lij or uij during each stage of Gaussian elimination (sometimes referred to as a partial sums accumulation). This advantage is based on the fact the product of two single precision floating point numbers is always computed with double precision arithmetic (at least in the C programming language). Because of this, the product aipapj suffers no loss of precision. If the product is accumulated in a double precision variable , there is no loss of precision during the entire inner product calculation. Therefore, one double precision variable can preserve the numerical integrity of the inner product.
Recalling the partial sum accumulation loop of the elimination-based procedure,
You will observe that truncation to single precision must occur each time a is updated unless both A and w are stored in double precision arrays.