4 LU Theorem

If you put minus the multipliers in a Gaussian elimination without partial pivoting in a lower triangular matrix $L$, with values 1 on the main diagonal, then

$$\fbox{$\displaystyle L U = A $}$$ (1)

That is the $LU$ theorem.

After $L$ and $U$ have been found through Gaussian elimination, you can solve $A\vec x=\vec b$ for any given right hand side vector $\vec b$ as follows: The system to solve is according to the theorem $LU\vec x=\vec b$. Temporarily call $U\vec x=\vec y$. Solve $L\vec y=\vec b$ using simple forward substitution. With $\vec y$ now known, $U\vec x=\vec y$ can be solved using back substitution as before:

$$\fbox{$\displaystyle L \vec y = \vec b \qquad U \vec x = \vec y $}$$ (2)

In case partial pivoting was needed, remember the row order interchanges you did and do the same interchanges with the coefficients of $\vec b$ before solving as above. You also need to do the row exchanges on the forming matrix $L$ while doing the GE. $LU$ will then be the matrix $A^{pp}$, equal to matrix $A$ except with its rows in the final order produced by all the row interchanges.

For large $n\times n$ matrices the number of computer operations (defined as 1 multiplication and 1 addition) to find $L$ and $U$ is approximately $\frac13n^3$. After that, to find a $\vec x$ given a $\vec b$ takes only about $n^2$ operations.

Some warnings. Normally speaking:

1. Never ever use Cramer's rule for anything but the tiniest of matrices. Multiplying out a determinant takes about $(n+1)!$ operations, which is gigantically larger than $\frac13n^3$ for everything except the tiniest $n$. Not to mention the possible round-off error growth with so many operations. And the risk of overflow and underflow (i.e. numbers getting outside the range that they can be stored on a computer.)
2. If you find the inverse matrix $A^{-1}$, you can simply find any $\vec x$ as $A^{-1}\vec b$. But that is normally a bad idea. One reason is that it takes $n^3$ operations to find $A^{-1}$, three times more than to find $L$ and $U$. And to evaluate $A^{-1}\vec b$ still takes $n^2$ operations. Also the additional operations tend to increase round-off error, [1].
3. A band matrix is a matrix in which the nonzero elements are restricted to a relatively narrow band around the main diagonal. Never ever use an LU decomposition library routine designed for a full matrix to solve a system with a band matrix. The waste in storage and computer effort to store all these zeros, and compute with them, would be gigantic. Use an LU decomposition subroutine designed for a band matrix.
4. Never ever find the inverse of a band matrix if you can help it. The inverse will be a full matrix. However, $L$ and $U$ will still be band matrices. (Partial pivoting may increase $U$ a bit to outside the band of $A$, by up to the width of $L$.)