clear
echo on
format compact
%This example demonstrates the PLU
%factorization of a square matrix A.
%This procedure is the Gaussian elimination
%process with row pivoting.

%First, we define our matrix A:
a=[1 3 5 7;1 2 3 4;5 3 2 1;3 4 2 6]
pause
%We determine the size of the matrix:
[n,m]=size(a)
%so we have n rows to operate on!  We
%will do the Gaussian elimination process
%a total of n-1 times.  We begin by initializing
%the p, l, and u matrices:
p=[1:n];
l=zeros(n,m);
u=a;

%Note that normally, to save space in memory,
%we would not save p, l and u separately, but
%rather l and u would overwrite the matrix a
%and the permutation matrix p would be saved
%as a single vector rather than a matrix.  We
%will not do this here, however, to make the
%operations clearer.
pause
%Now let's begin:
for i=1:n-1
 %First we check to be sure that the largest
 %element in the column to be operated on (the
 %i'th column) is in the ith row:
 [junk,ibig]=max(abs(u(i:n,i)));
 ibig=ibig+(i-1)
 %so we trade the two relevant rows of u and
 %update the pivot matrix p:
 temp=u(i,i:m);
 u(i,i:m)=u(ibig,i:m);
 u(ibig,i:m)=temp;
 temp=p(i);
 p(i)=p(ibig);
 p(ibig)=temp;
 pause
 %Let's look at the matrices u and p:
 u
 pmat=zeros(m,m);
 for k=1:m
  pmat(p(k),k)=1;
 end
 pmat
 %Note that we now have the correct row in our
 %matrix u, and that the permutation matrix p
 %has also been adjusted.
 pause
 %At the same time that we switch rows in u, we
 %also have to switch rows in l:
 temp=l(i,:);
 l(i,:)=l(ibig,:);
 l(ibig,:)=temp;
 %Now for the Gaussian elimination part.  First
 %we will determine the multipliers:
 mult=u(i+1:n,i)./u(i,i);
 %This array of multipliers becomes the i'th
 %column of the lower triangular matrix l:
 l(i+1:n,i)=mult;
 %and now we do the reduction:
 for k=i+1:n,
  u(k,:)=u(k,:)-mult(k-i)*u(i,:);
 end
 %Our matrix u now looks like:
 u
 %and we are ready for the next row!
 pause
end
%To finish, we add a diagonal of ones to the
%matrix l:
l=l+eye(n)
%and the other two matrices are:
pmat
u
pause
%We may verify the decomposition by comparing
%to our original matrix a:
pmat*l*u
%which matches a:
a
%And thus we can use p, l and u to solve for
%the solution to any right hand side vector b.
pause
%As a final note, the matlab command [L,U,P]=lu(a)
%does the same thing as this program, except that
%the matrix P is the inverse (transpose for
%this kind of permutation matrix) of the matrix
%pmat given in this program With the matlab
%decomposition, we get that a = P'*L*U instead
%of p*l*u.
[L,U,P]=lu(a)
A=P'*L*U