function [L U] = ilub(A)
%ilub incomplete LU factorization
%   [L U] = ilub(A) returns factors that may be useful as
%   preconditioners for a system of linear equations being solved by
%   iterative methods such as GMRES.

[m n] = size(A);
if m ~= n
   error('Matrix must be square.');
end

% store L and U in A
for i = 2:n
   for j = 1:i-1
      if A(i,j) ~= 0.0
         % store nonzero multiplier
         A(i,j) = A(i,j)/A(j,j);
         for k = j+1:n
            if A(i,k) ~= 0.0
               A(i,k) = A(i,k) - A(i,j)*A(j,k);
            end
         end
      end
   end
end

U = triu(A);
L =  tril(A,-1);
% make the unit diagonal
for i = 1:n
	L(i,i) = 1;
end;

