MatrixLibrary | A matrix library written in C | Math library

When you perform a Gaussian elimination, you swap rows and repeatedly subtract a multiple of one row from another to produce an upper triangular form.

When you do this on a system of equations or an "augmented matrix", you do not use any information from the result column. The decisions about which rows to swap and which to subtract with what multiplier would be exactly the same no matter what numbers are in the result column.

Because the "result column" is not used, you can perform the same procedure on a normal square matrix. Since the operations don't change the determinant (if you negate one row whenever you swap), you end up with an upper triangular matrix with the same det as the original.

The MIT author calls a function lu to do this in the example near the start. This does L-U decomposition on the matrix, which returns the Gaussian-eliminated matrix in the 'U' part: https://en.wikipedia.org/wiki/LU_decomposition.

L-U decomposition is pretty cool. It's like doing Gaussian elimination to solve all systems with the same "matrix part" all at once, which again you can do because the process doesn't need to see the result column at all.

Starting with a matrix M, you get L and U such that LU = M. That means, if you want to solve:

Mx = y

... where (x an y are column vectors), you have:

LUx = y

Solve Lv=y, which is easy (just substitution) because L is lower-triangular. Then you have:

Ux = v

... which is easy to solve because U is upper-triangular.