4.1 Gaussian Elimination

Approaches to LU decomposition which systematically capture the intermediate results of Gaussian elimination often differ in the order in which A is forced into upper triangular form. The most common alternatives are to eliminate the subdiagonal parts of A either one row at a time or one column at a time. The calculations required to zero a complete row or a complete column are referred to as one stage of the elimination process.

The effects of the kth stage of Gaussian elimination on the A matrix are summarized by the following equation.

$$a_{ij}^{(k+1)}=a_{ij}^{\left(k\right)}-\left(\frac{a_{ik}^{\left(k\right)}}{a_{kk}^{\left(k\right)}}\right)a_{ij}^{\left(k\right)},\mbox{ where }i,j \gt k$$ (43)

The notation $a_{ij}^{(k)}$ means the value of $a_{ij}$ produced during the kth stage of the elimination procedure. In Equation 43, the term $\frac{a_{ik}^{\left(k\right)}}{a_{kk}^{\left(k\right)}}$ (sometimes referred to as a multiplier) captures the crux of the elimination process. It describes the effect of eliminating element $a_{ik}$ on the other entries in row i during the kth stage of the elimination. In fact, these multipliers are the elements of the lower triangular matrix L, i.e.

$$l_{ik}=\frac{a_{ik}^{\left(k\right)}}{a_{kk}^{\left(k\right)}}$$ (44)

Algorithm 1 implements Equation 43 and Equation 44 and computes the LU decomposition of an m × n matrix A. It is based on Algorithm 4.2–1 of  Golub and Van Loan[2] .

Algorithm 1: Gaussian Elimination
for $k=1,\cdots,\min(m-1,n)$
   for $j=k+1,\cdots,n$
   for $i=k+1,\cdots,m$
      for $j=k+1,\cdots,n$
         $a_{ij}=a_{ij}-\alpha w_{j}$

The algorithm overwrites aij with lij when i < j. Otherwise, aij is overwritten by uij. The algorithm creates a matrix U that is upper triangular and a matrix L that is unit lower triangular. Note that a working vector w of length n is required by the algorithm.