%% Using gmres Preconditioner with Restart
% Using a preconditioner with restarted gmres.
% Load west0479, a real 479-by-479 nonsymmetric sparse matrix.

load west0479
A = west0479;

%
% Define b so that the true solution is a vector of all ones.

b = sum(A,2);

%
% Construct an incomplete LU preconditioner with a drop tolerance of 1e-6.

[L,U] = ilu(A,struct('type','ilutp','droptol',1e-6));

% The benefit to using restarted gmres is to limit the amount of memory required to execute the method. Without restart, gmres requires maxit vectors of storage to keep the basis of the Krylov subspace. Also, gmres must orthogonalize against all of the previous vectors at each step. Restarting limits the amount of workspace used and the amount of work done per outer iteration.
%
% Execute gmres(3), gmres(4), and gmres(5) using the incomplete LU factors as preconditioners. Use a tolerance of 1e-12 and a maximum of 20 outer iterations.

tol = 1e-12;
maxit = 20;
[x3,fl3,rr3,it3,rv3] = gmres(A,b,3,tol,maxit,L,U);
[x4,fl4,rr4,it4,rv4] = gmres(A,b,4,tol,maxit,L,U);
[x5,fl5,rr5,it5,rv5] = gmres(A,b,5,tol,maxit,L,U);
fl3
fl4
fl5

% fl3, fl4, and fl5 are all 0 because in each case restarted gmres drives the relative residual to less than the prescribed tolerance of 1e-12.
%
% The following plot shows the convergence history of each restarted gmres method. gmres(3) converges at outer iteration 5, inner iteration 3 (it3 = [5, 3]) which would be the same as outer iteration 6, inner iteration 0, hence the marking of 6 on the final tick mark.

semilogy(1:1/3:6,rv3/norm(U\(L\b)),'-o');
h1 = gca;
h1.XTick = (1:6);
title('gmres(N) for N = 3, 4, 5')
xlabel('Outer iteration number');
ylabel('Relative residual');
hold on
semilogy(1:1/4:3,rv4/norm(U\(L\b)),'-o');
semilogy(1:1/5:2.8,rv5/norm(U\(L\b)),'-o');
yline(tol,'r--');
hold off
legend('gmres(3)','gmres(4)','gmres(5)','Tolerance')
grid on

% In general the larger the number of inner iterations, the more work gmres does per outer iteration and the faster it can converge.
% Copyright 2019 The MathWorks, Inc.