Now we can address the general LU-solver. Assuming the lower and
upper triangular factorization, we will solve
. We will assume here that L has 1s on the
diagonal; call the other non-zero entries
of L 
. The entries of U will be 
.
Arithmetic tells us (Crout's algorithm)
Convince yourself that the terms appearing on the right-hand-side
of these equations are available by the time they are needed.
You will also see that each 
is used only once; the
can be stored in their place in A.
Partial pivoting can be implemented (i.e., row interchange). This
means we need to keep track of any permutation that is performed.
This is implemented via the index array in the code.
The algorithm below also scales each row by the largest entry in that row
- this is the implicit scaling we saw to be effective in the notes.
The algorithm below first factors A into L and U (ludcmp). Given this decomposition, lubksb solves the system Ax=b. Note that if you need to solve several systems of the form Ax=b with different b's, you need call the decomposition only once, and then make repeated calls to the back substitution routine.
A quick estimate shows the decomposition takes about
operations; the total count for solving a system is 
.
If you want to invert A, find the inverse column-by-column. That
is, solve N times with the right-hand-side being the unit vector
- 0s, with a 1 in the kth entry.