function [l,u]=gauss_elim(a);
% Usage: [l,u]=gauss_elim(a);
% Reduces matrix a to upper triangular form using Gaussian elimination
% Input:
% a = input matrix
% Output:
% l = lower triangular transformation of a
% u = upper triangular transformation of a

[n,m]=size(a);
for k=1:n-1 % Loop over columns
    if a(k,k)==0
       break; % Stop if pivot is zero
    end
    a(k+1:n,k)=a(k+1:n,k)/a(k,k); % Compute multipliers for column k
    for j = k+1:n
	a(k+1:n,j) = a(k+1:n,j) - a(k+1:n,k)*a(k,j);		
    end
end
u=triu(a); % Upper trianguler matrix
l=tril(a,-1)+eye(n); % Unit lower triangular matrix