LU decomposition

% Yuri Lvov code

% direct implementation of LU factorization from page 69

% of the book by Heath;

% Thu Sep 19 14:45:19 EDT 2002

N=5; %Size of MAtrix

disp('Compare buid in and the one I wrote based on direct implementation');

A=rand(N,N); SaveA=A;

m=eye(N,N);

U=zeros(N,N);

disp('matrix A'); disp(A);

tic;

[l,u]=lu(A);

toc

tic;

for k=1 :N-1, %Loop over columns

if A(k,k)==0 disp('Zero pivot, sorry '); return; end;

for i=k+1:N, %Compute multipies

m(i,k)=A(i,k)/A(k,k); % I can store result in A,

% here m is for psyco convenience

end;

for j=k+1:N,

for i = k+1:N,

A(i,j)=A(i,j)- m(i,k)*A(k,j);

end;

end;

end;

toc;

disp('resulting A');disp(A);

disp('resulting m');disp(m);

for i=1:N,

for j=i:N,

U(i,j)=A(i,j); end;

end;

disp(' Check Matlab accuracy'); disp(norm(SaveA-l * u ));

disp(' Check Heath accuracy'); disp(norm(SaveA-m * U ));