  function [L,U] = GE(A);
%
% Pre:
%   A      n-by-n
%
% Post:
%   L      n-by-n unit lower triangular with |L(i,j)|<=1.
%   U      n-by-n upper triangular.
%          A = LU  
 
    [n,n] = size(A); 
    for k=1:n-1
       A(k+1:n,k) = A(k+1:n,k)/A(k,k);
       A(k+1:n,k+1:n) = A(k+1:n,k+1:n) - A(k+1:n,k)*A(k,k+1:n);
    end
	L = eye(n,n) + tril(A,-1);
    U = triu(A);

         
