function [L,U] = gelu(A)
% LU-factorization based on Gaussian elimination
% No pivoting: UNSTABLE!
[m,n] = size(A); L = eye(n);
% Forward elimination
for i=1:n-1
  if (A(i,i) == 0.0), error('Zero pivot'); end
  for j=i+1:n
    fac = A(j,i)/A(i,i); L(j,i) = fac;
    for k=i+1:n
      A(j,k) = A(j,k) - fac*A(i,k); 
    end
  end
end
U = triu(A);