%% Using gmres Preconditioner without Restart
% Examine the effect of using a preconditioner matrix with non-restarted gmres to solve a linear system.
% Load west0479, a real 479-by-479 nonsymmetric sparse matrix.

load west0479
A = west0479;

%
% Define b so that the true solution to \textrm{Ax}=b is a vector of all ones.

b = sum(A,2);

%
% Set the tolerance and maximum number of iterations.

tol = 1e-12;
maxit = 20;

%
% Use gmres to find a solution at the requested tolerance and number of iterations. Specify five outputs to return information about the solution process:
% x is the computed solution to A*x = b.
% fl0 is a flag indicating whether the algorithm converged.
% rr0 is the relative residual of the computed answer x.
% it0 is a two-element vector [outer inner] indicating the inner and outer iteration numbers when x was computed.
% rv0 is a vector of the residual history for \left\|b-\textrm{Ax}\right\|.

[x,fl0,rr0,it0,rv0] = gmres(A,b,[],tol,maxit);
fl0
rr0
it0

% fl0 is 1 because gmres does not converge to the requested tolerance 1e-12 within the requested 20 iterations. The best approximate solution that gmres returns is the last one (as indicated by it0(2) = 20). MATLAB stores the residual history in rv0.
%
% To aid with the slow convergence, you can specify a preconditioner matrix. Since A is nonsymmetric, use ilu to generate the preconditioner M=L\;U. Specify a drop tolerance to ignore nondiagonal entries with values smaller than 1e-6. Solve the preconditioned system M^{-1} A\;x=M^{-1} b by specifying L and U as inputs to gmres. Note that when you specify a preconditioner, gmres calculates the residual norm of the preconditioned system for the outputs rr1 and rv1.

[L,U] = ilu(A,struct('type','ilutp','droptol',1e-6));
[x1,fl1,rr1,it1,rv1] = gmres(A,b,[],tol,maxit,L,U);
fl1
rr1
it1

% The use of an ilu preconditioner produces a relative residual less than the prescribed tolerance of 1e-12 at the sixth iteration. The first residual rv1(1) is norm(U\(L\b)), where M = L*U. The last residual rv1(end) is norm(U\(L\(b-A*x1))).
%
% You can follow the progress of gmres by plotting the relative residuals at each iteration. Plot the residual history of each solution with a line for the specified tolerance.

semilogy(0:length(rv0)-1,rv0/norm(b),'-o')
hold on
semilogy(0:length(rv1)-1,rv1/norm(U\(L\b)),'-o')
yline(tol,'r--');
legend('No preconditioner','ILU preconditioner','Tolerance','Location','East')
xlabel('Iteration number')
ylabel('Relative residual')

% Copyright 2020 The MathWorks, Inc.