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] .
| for $k=1,\cdots,\min(m-1,n)$ |
| for $j=k+1,\cdots,n$ |
| $w_{j}=a_{kj}$ |
| for $i=k+1,\cdots,m$ |
| $\alpha=\displaystyle\frac{a_{ik}}{a_{kk}}$ |
| $a_{ik}=\alpha$ |
| 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.