%% Using bicgstabl with Preconditioner
% Examine the effect of using a preconditioner matrix with bicgstabl 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 bicgstabl to find a solution at the requested tolerance and number of iterations. Specify five 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\|.

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

% fl0 is 1 because bicgstabl does not converge to the requested tolerance 1e-12 within the requested 20 iterations. In fact, the behavior of bicgstabl is so poor that the initial guess x0 = zeros(size(A,2),1) is the best solution and is returned, as indicated by it0 = 0.
%
% 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 {A\;M}^{-1} M\;x=b by specifying L and U as inputs to bicgstabl.

setup = struct('type','ilutp','droptol',1e-6);
[L,U] = ilu(A,setup);
[x1,fl1,rr1,it1,rv1] = bicgstabl(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 second iteration. The output rv1(1) is norm(b), and the output rv1(end) is norm(b-A*x1).
%
% You can follow the progress of bicgstabl 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 2015 The MathWorks, Inc.
%
%