function [l,u] = simple_lu ( a )
%
%  [l,u] = simple_lu(a) applies the LU factorization 
%    to the matrix in a, without pivoting.
%
%  example: a = randn(3,3)
%           [l,u] = simple_lu(a)
%           l*u - a
%
n = size(a,1);
b = a;
for j=1:n-1
  for i=j+1:n
    b(i,j) = b(i,j)/b(j,j);
  end;
  for i=j+1:n
    for k=j+1:n 
      b(i,k) = b(i,k) - b(i,j)*b(j,k);
    end;
  end;
end;
% the lower part of the matrix b
% contains l, the upper part contains u.
l = diag(ones(1,n)) + tril(b,-1);
u = triu(b);