%% Using pcg with Preconditioner
% Examine the effect of using a preconditioner matrix with pcg to solve a linear system.
%
% Create a symmetric positive definite, banded coefficient matrix.

A = delsq(numgrid('S',102));

%
% Define b for the right-hand side of the linear equation.

b = ones(size(A,1),1);

% Set the tolerance and maximum number of iterations.

tol = 1e-8;
maxit = 100;

%
% Use pcg 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] = pcg(A,b,tol,maxit);
fl0
rr0
it0

% fl0 is 1 because pcg does not converge to the requested tolerance of 1e-8 within the requested 100 iterations.
%
% To aid with the slow convergence, you can specify a preconditioner matrix. Since A is symmetric, use ichol to generate the preconditioner M=L\;L^T. Solve the preconditioned system by specifying L and L' as inputs to pcg.

L = ichol(A);
[x1,fl1,rr1,it1,rv1] = pcg(A,b,tol,maxit,L,L');
fl1
rr1
it1

% The use of an ichol preconditioner produces a relative residual less than the prescribed tolerance of 1e-8 at the 79th iteration. The output rv1(1) is norm(b) and rv1(end) is norm(b-A*x1).
%
% Now, use the michol option to create a modified incomplete Cholesky preconditioner.

L = ichol(A,struct('michol','on'));
[x2,fl2,rr2,it2,rv2] = pcg(A,b,tol,maxit,L,L');
fl2
rr2
it2

% This preconditioner is better than the one produced by the incomplete Cholesky factorization with zero fill for the coefficient matrix in this example, so pcg is able to converge even quicker.
%
% You can see how the preconditioners affect the rate of convergence of pcg by plotting each of the residual histories starting from the initial estimate (iterate number 0). Add 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')
semilogy(0:length(rv2)-1,rv2/norm(b),'-o')
yline(tol,'r--');
legend('No Preconditioner','Default ICHOL','Modified ICHOL','Tolerance','Location','East')
xlabel('Iteration number')
ylabel('Relative residual')

% Copyright 2019 The MathWorks, Inc.