function [L, U] = ilup(A,p) 
n = size(A,1) ;
L = speye(n,n); 
U = sparse(n,n);
levels = sparse(n,n); 
for i = 1:n
  w = A(i,:);
  jr = find(w) ;
  levi = sparse(1,n) ; 
  levi(jr) = 1;
  j = 0;
  while (j < i)
    % select smallest column index among jw(k), k=j+1, ..., lenl. 
    j = j+min(find(w(j+1:n))) ;
    if (j >= i)
      break;
    end 
    jlev = levi(j) ;
    if (jlev <= (p-1)) 
      %      combine current row and  row j.  
      [ir,jr,v] = find(U(j,:)); 
      fact = w(j)/v(1); 
      for k=1:length(jr) 
        jj = jr(k) ;
        temp = levels(j,jj) + jlev ; 
        if (levi(jj) == 0 ) 
          if (temp <= (p-1)) 
            w(jj) =  - fact*v(k); 
            levi(jj) = temp;
          end
        else 
          w(jj) = w(jj) - fact*v(k); 
          levi(jj) = min(levi(jj),temp); 
        end
      end 
      w(j) = fact;   
    end
  end
  [iw, jw, s] = find(w) ;
  [il,jl,ll] = find(levi); 
  for k=1:length(jw)  
    j = jw(k) ;
    if (j < i)
      L(i,j) = s(k)  ; 
    else 
      U(i,j) = s(k) ;
      levels(i,j) = ll(k) ; 
    end
  end   
end