function [ L,U ] = LUdecomp(A)
%LUdecomp: decompose square matrix A as A=LU
% where L is lower triag and U is upper triag
    [m, n]=size(A);
    if m ~= n, error('Input must be a square matrix.'), end
    L=zeros(n); U=zeros(n);

    % remember A^(0) is A
    AofK = A;  
    for k = 1:n
        % at this point AofK is A^(k-1)
        for j = k:n
            U(k,j) = AofK(k,j);
        end
        % check that we don't divide by zero(!)
        if U(k,k) == 0
            error('** A^(k-1)_{k,k}==0 in LU decomp')
        end
        for i = k:n
           L(i,k) = AofK(i,k)/U(k,k);
        end
        % now modify AofK so that we can use it in the
        % next iteration
        for i = k:n
            for j = k:n
                AofK(i,j) = AofK(i,j) - L(i,k)*U(k,j);
            end
        end
    end

end
