function [L,U]=ILU0(A)
% ILU0 incomplete ILU factorization with zero fill-in
%   [L,U]=ILU0(A); computes an ILU factorization of the matrix
%   A with zero fill-in.

n=size(A,1);
for k=1:n-1
  for i=k+1:n
    if A(i,k)~=0
      A(i,k)=A(i,k)/A(k,k);
      for j=k+1:n
        if A(i,j)~=0
          A(i,j)=A(i,j)-A(i,k)*A(k,j);
        end
      end
    end
  end
end
L=speye(size(A))+tril(A,-1);
U=triu(A,0);
