function [L,A] = LU_factor(A,n)
% LU factorization of an n by n matrix A
% using Gauss elimination without pivoting.
% A is factored as A = L*U
% Output:
% L is lower triangular with the main diagonal part = 1s.
% U is upper triangular & is stored in the original mtx A.
% Its lower triangular part has not been zeroed out to get U.
% We don't bother to do that here since back-substitution
% will not even use the lower triangular part of U.
% K. Ming Leung, 01/26/03, 
% K. Ming Leung, 02/21/04 (version 2) No unnecessary computation &
% more efficient use of indexing.

L = eye(n);
for k = 1:n-1
    if (A(k,k) == 0) Error('Pivoting is needed!'); end
    L(k+1:n,k) = A(k+1:n,k)/A(k,k);
    i = k+1:n;
    A(i,i) = A(i,i) - L(i,k)*A(k,i);
end
