%% 
% Here is a 1000x1000 matrix of density around 0.5%.
rng(9)
A = 0.6*speye(1000) + sprand(1000, 1000, 0.005, 1/10000);

%%
% Without a preconditioner, GMRES takes a large number of iterations. 
b = rand(1000, 1); 
[x, ~, ~, ~, resid_plain] = gmres(A, b, 50, 1e-10, 6);

figure(1)
semilogy(resid_plain,'-')
grid on
xlabel('iteration number')
ylabel('residual norm')
title('Unpreconditioned GMRES')

%%
% This version of incomplete LU factorization simply prohibits fill-in
% for the factors, freezing the sparsity pattern of the approximate
% factors.
[L,U] = ilu(A);

figure(2)
subplot(121)
spy(A)
title('A')
subplot(122)
spy(L)
title('L')

%%
% It does _not_ produce a true factorization of A.
norm(full(A - L*U))  

%%
% The actual preconditioning matrix is M = LU. However, the
% |gmres| function allows setting the preconditioner by giving the factors
% independently.
[x, ~, ~, ~, resid_prec] = gmres(A, b, [], 1e-10, 300, L, U);

%%
% The preconditioning is fairly successful in this case.
figure(3)
semilogy(resid_prec,'-')
grid on
xlabel('iteration number')
ylabel('residual norm')
title('Precondtioned GMRES')