LU decomposition (with pivoting)

In a first course you ran Gaussian elimination to solve Ax = b: you subtracted multiples of one row from the rows below until the matrix was upper triangular. The advanced insight is that elimination is itself a factorization. Every row operation you did was a multiplication by a simple lower-triangular matrix, and bookkeeping those operations gives A = L U, where L is unit lower triangular (ones on the diagonal, the elimination multipliers below) and U is the upper triangular result.

Plain elimination breaks the moment a pivot is zero, and it behaves badly when a pivot is merely tiny, because dividing by a small number magnifies rounding error. The cure is partial pivoting: before eliminating in a column, swap rows to bring the largest-magnitude entry to the pivot. Recording the swaps in a permutation matrix P gives the practical statement P A = L U. Now every multiplier in L has magnitude at most 1.

This is what production solvers actually compute. Factor once at cost about (2/3) n^3 flops, then solve A x = b for any right-hand side by a fast forward solve L y = P b and back solve U x = y, each costing only n^2. The same factorization yields det(A) = (+/-1) times the product of the diagonal of U, the sign coming from the parity of the row swaps.

Honesty caveat: partial pivoting is stable in practice for nearly all matrices, but it is not unconditionally stable. There exist contrived matrices (the classic example grows like 2^n) where the entries of U blow up. These almost never appear in real problems, which is why partial pivoting remains the default rather than the more expensive complete pivoting.