 function [L,U] = ilu0(A) 
%---------------------------------------------
% function [L,U] = ilu0 (A) 
% ILU factorization of A. Uses ikj variant of GE
%---------------------------------------------
n = size(A,1) ;
for i=2:n 
    [ii,jj,rr] = find(A(i,:));
    p = 1; 
    k = jj(p) ;
    while (k < i)    
        piv = A(i,k)/A(k,k); 
        A(i,k) = piv;
        p=p+1;
        u = jj(p:end);
%%QUESTION: this will be slow. Why? Alternatives? 
        A(i,u) = A(i,u)-piv*A(k,u);
        k = jj(p) ;
     end  %% while
 end 
 L = tril(A,-1) + spdiags(ones(n,1),0,n,n) ;
 U = triu(A);
