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 k^{th} 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)},\operatorname{where}i,j>k$ | (27) |

The notation a${}_{ij}^{(k)}$ means the value of a_{ij} produced during the
k^{th} stage of the elimination procedure. In Equation
27, 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
k^{th} 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)}}$ | (28) |

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

for $k=1,\cdots,\operatorname{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 a_{ij} with l_{ij}
when i < j. Otherwise, a_{ij} is overwritten by
u_{ij}. 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.