clear
echo on
%In this example we demonstrate the use
%of PLU decomposition of a matrix 'a' to
%solve a system of equations.  First we
%set up the matrix.  We will use a 7x7
%hilbert matrix (a relatively ill-conditioned
%matrix).

n=7;
a=zeros(n,n);
for i=1:n
 for j=1:n
  a(i,j)=1/(i+j-1);
 end;
end;
pause
%Now let's display this matrix:
a
pause
%We also want to set up the right hand side:
x=[1:n]';
b=a*x
%Which gives us the rhs which has as its solution
%the array x = [1:n]'.
pause

%Now let's solve the system using LU decomposition:
[l u pt]=lu(a);
%which generates the lower triangular and upper
%triangular matrices l and u, and the transpose of
%the permutation matrix pt=p'.  The original matrix
%a can be recovered by a = p * l * u.  For more
%information, try typing in 'help lu'.  Let's
%examine each of these.
pause
%The permutation matrix is given by:
p=pt'
%Note that it is only a permuted identity matrix - all
%it does is trade rows of the matrix a, corresponding
%to pivoting.
pause
%The matrix p has the property that it is orthogonal,
%e.g. that p'*p = eye.  We can show this:
p'*p
pause
%The matrix l is a lower triangular matrix with ones on
%the diagonal.  The elements are simply the multipliers
%in a gaussian elimination process:
l
pause
%The matrix u is the upper triangular matrix - what is
%left over after gaussian elimination is complete:
u
pause
%We can use the matrices to solve the original system
%for arbitrary rhs b by solving the series of problems:
%  p * z = b
%  l * y = z
%  u * x = y
%This is much faster than doing gaussian elimination
%over from scratch for every b.  Let's show this:
pause
flops(0)
z=p'*b;
y=l\z;
x=u\y;
%So the answer is x:
x
%and the total number of operations in the back and
%forward substitution was:
subflops=flops
pause
%We can compare this to solving for x directly for
%each rhs b:
flops(0)
x=a\b
elimflops=flops
%So the answer is the same, but the number of operations
%is higher.  The discrepancy would be much larger if
%n were larger.
echo off