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

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

%
% 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 = 100;

%
% Use symmlq 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*x = 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\|.
% rvcg0 is a vector of the conjugate gradient residual history for \left\|A^T A\;x-A^T b\right\|.

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

% fl0 is 1 because symmlq does not converge to the requested tolerance 1e-12 within the requested 100 iterations.
%
% To aid with convergence, you can specify a preconditioner matrix. Since A is symmetric, use ichol to generate the preconditioner M=L\;L^T. Specify the 'ict' option to use incomplete Cholesky factorization with threshold dropping, and specify a diagonal shift value of 1e-6 to avoid nonpositive pivots. Solve the preconditioned system by specifying L and L' as inputs to symmlq.

setup = struct('type','ict','diagcomp',1e-6,'droptol',1e-14);
L = ichol(A,setup);
[x1,fl1,rr1,it1,rv1,rvcg1] = symmlq(A,b,tol,maxit,L,L');
fl1
rr1
it1

% The use of an ichol preconditioner greatly improves the numerical properties of the problem, and symmlq is able to converge quickly. The output rv1(1) is norm(b), and the output rv1(end) is norm(b-A*x1).
%
% You can follow the progress of symmlq by plotting the relative residuals at each iteration. Plot the conjugate gradient residual history of each solution with a line for the specified tolerance.

semilogy(0:length(rvcg0)-1,rvcg0/norm(b),'-o')
hold on
semilogy(0:length(rvcg1)-1,rvcg1/norm(b),'-o')
yline(tol,'r--');
legend('No preconditioner','ICHOL preconditioner','Tolerance','Location','SouthEast')
xlabel('Iteration number')
ylabel('Relative residual')

% Copyright 2019 The MathWorks, Inc.