function [L,A] = lunp(A)
% [L,U] = lunp(A)
% Gaussian elimination, no pivoting: find L, U such that A=L*U
% where L is lower triangular with 1 on diagonal,
%       U is upper triangular, nonzero elements on diagonal
%
% U = lunp(A) only returns the matrix U as defined above  
m = size(A,1); n = size(A,2); 
L = eye(m);
if isa(A,'sym')
  L = sym(L);
end
for k=1:min(m,n)
  if k<m
    if A(k,k)==0
      warning('LU-decomposition does not exist');
      L=[]; A=[];
      return;
    end  
    L(k+1:end,k) = A(k+1:end,k)/A(k,k);    % compute multipliers
    A(k+1:end,k) = zeros(m-k,1);
    if k<n
      A(k+1:end,k+1:end) = A(k+1:end,k+1:end) - L(k+1:end,k)*A(k,k+1:end);
    end
  end
end

if nargout<2
  L=A;
end
