function [ x,r,S,T,iter ] = ilu_iter_solv( A,x0,b,tol,droptol )

setup.type = 'ilutp';
setup.droptol = droptol;
[L,U] = ilu(sparse(A),setup);
U = full(U);
L = full(L);

S = L*U; 
T = S-A;
r = b-A*x0;

iter = 0;
while norm(r) > tol
iter = iter + 1;    
x1 = S\(T*x0+b);
r = b-A*x1;
x0 = x1;
end

x = x1;

end

