function [L U]= luf(A)
% LU factorization of a square matrix: A= LU

[m n]=size(A);

% initialize L and U
L=zeros(n); U=zeros(n);

% start the factorization
for k=1:n  
    % compute L(k,k) and U(k,k)
    p= dot( L(k,1:k-1), U(1:k-1,k)' );
    L(k,k)= 1.0;
    U(k,k)= ( A(k,k)-p )/L(k,k);
    
    % determine the kth column of L
    q= L(k+1:n,1:k-1)*U(1:k-1,k); 
    L(k+1:n, k) = (A(k+1:n, k) - q)/U(k,k);
   
    % calculate the kth row of U
    r = L(k,1:k-1)*U(1:k-1,k+1:n);
    U(k, k+1:n)= (A(k, k+1:n) - r)/L(k,k);
end