Matrices are used widely in Numerical Computing, many algorithms
involve decomposing a matrix, A into an equivalent pair L and
U, where L and U are lower and upper triangular
matrices respectively and .
Consider the following system of equations,
here
so
and
Both L and U can be stored in the same matrix; since the diagonal of L is always unity it does not have to be store, the combined LU matrix can be stored as
Now, given that
then we solve for where
(recall that the diagonal of L will be all 1s) and then solve
The advantage of this method is that the triangular system of equations
is very easy to solve.
can be solved by forward
substitution and
by
backward substitution.
is calculated by forward
substitution as follows :
Once is known
can easily be obtained by
backwards substitution:
The LU decomposition algorithm is as follows:
For each row k from 1 to n-1 in order,
A(k,k)
are divided by the diagonal element.A(i,j)
in the submatrix of A
below and to the right of the diagonal
element A(k,k)
, ie, A(k+1:n,k+1:n)
is modified by
subtracting A(i,k)*A(k,j)
The algorithm will not work if the A
matrix is singular, that is,
if any of the rows of A
are all zero.
This algorithm omits the issue of pivoting so can be unstable - your program will work for the 3 by 3 matrix given here.
The an HPF program with subroutines to perform LU decomposition and forward and backwards substitution.