%% Using lsqr with Preconditioner
% Examine the effect of using a preconditioner matrix with lsqr 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 is a vector of all ones.

b = sum(A,2);

%
% Set the tolerance and maximum number of iterations.

tol = 1e-12;
maxit = 20;

%
% Use lsqr to find a solution at the requested tolerance and number of iterations. Specify six outputs to return information about the solution process:
% x0 is the computed solution to A*x0 = b.
% fl0 is a flag indicating whether the algorithm converged.
% rr0 is the residual of the computed answer x0.
% it0 is the iteration number when x0 was computed.
% rv0 is a vector of the residual history for \left\|\textrm{Ax}-b\right\|.
% lsrv0 is a vector of the least squares residual history.

[x0,fl0,rr0,it0,rv0,lsrv0] = lsqr(A,b,tol,maxit);
fl0
rr0
it0

% Since fl0 = 1, the algorithm did not converge to the specified tolerance within the maximum number of iterations.
%
% 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 \tilde{A} y={\textrm{AM}}^{-1} \left(M\;x\right)=b by specifying L and U as inputs to lsqr.

setup = struct('type','ilutp','droptol',1e-6);
[L,U] = ilu(A,setup);
[x1,fl1,rr1,it1,rv1,lsrv1] = lsqr(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 13th iteration. The output rv1(1) is norm(b), and the output rv1(end) is norm(b-A*x1).
%
% You can follow the progress of lsqr 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(b),'-o')
yline(tol,'r--');
legend('No preconditioner','ILU preconditioner','Tolerance','Location','East')
xlabel('Iteration number')
ylabel('Relative residual')

%
% Copyright 2019 The MathWorks, Inc.