% Gaussian elimination without pivoting
% function [l,u,err]=genp(a)
%
% inputs:
%   a - square matrix
%
% outputs:
%   l - unit lower triangular matrix
%   u - upper triangular matrix
%       (on output l*u = a, modulo roundoff)
%   err - array of diagnostic information:
%     err(1) = norm(l*u-a,1)/norm(a,1), a measure of backward error
%     err(2) = cond(a)
%     err(3) = cond(l)
%     err(4) = cond(u)
%
function [l,u,err]=genp(a)
n=max(size(a));
aa=a;
for i=1:n-1,
%  eliminate
   aa(i+1:n,i) = aa(i+1:n,i)/aa(i,i);
   aa(i+1:n,i+1:n) = aa(i+1:n,i+1:n) - aa(i+1:n,i)*aa(i,i+1:n);
end
l=eye(n);l=l+tril(aa,-1);
u=triu(aa,0);
err=[norm(l*u-a,1)/norm(a,1); cond(a); cond(l); cond(u)];
