%% Using ichol as a Preconditioner
% This example shows how to use an incomplete Cholesky factorization as a preconditioner to improve convergence.
% Create a symmetric positive definite matrix, A.

N = 100;
A = delsq(numgrid('S',N));

%
% Create an incomplete Cholesky factorization as a preconditioner for pcg. Use a constant vector as the right hand side. As a baseline, execute pcg without a preconditioner.

b = ones(size(A,1),1);
tol = 1e-6;
maxit = 100;
[x0,fl0,rr0,it0,rv0] = pcg(A,b,tol,maxit);

%
% Note that fl0 = 1 indicating that pcg did not drive the relative residual to the requested tolerance in the maximum allowed iterations. Try the no-fill incomplete Cholesky factorization as a preconditioner.

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

%
% fl1 = 0, indicating that pcg converged to the requested tolerance and did so in 59 iterations (the value of it1). Since this matrix is a discretized Laplacian, however, using modified incomplete Cholesky can create a better preconditioner. A modified incomplete Cholesky factorization constructs an approximate factorization that preserves the action of the operator on the constant vector. That is, norm(A*e-L*(L'*e)) will be approximately zero for e = ones(size(A,2),1) even though norm(A-L*L','fro')/norm(A,'fro') is not close to zero. It is not necessary to specify type for this syntax since nofill is the default, but it is good practice.

opts.type = 'nofill';
opts.michol = 'on';
L2 = ichol(A,opts);
e = ones(size(A,2),1);
norm(A*e-L2*(L2'*e))
[x2,fl2,rr2,it2,rv2] = pcg(A,b,tol,maxit,L2,L2');

%
% pcg converges (fl2 = 0) but in only 38 iterations. Plotting all three convergence histories shows the convergence.

semilogy(0:maxit,rv0./norm(b),'b.');
hold on
semilogy(0:it1,rv1./norm(b),'r.');
semilogy(0:it2,rv2./norm(b),'k.');
legend('No Preconditioner','IC(0)','MIC(0)');

%
% The plot shows that the modified incomplete Cholesky preconditioner creates a much faster convergence.
% You can also try incomplete Cholesky factorizations with threshold dropping. The following plot shows convergence of pcg with preconditioners constructed with various drop tolerances.

L3 = ichol(A, struct('type','ict','droptol',1e-1));
[x3,fl3,rr3,it3,rv3] = pcg(A,b,tol,maxit,L3,L3');
L4 = ichol(A, struct('type','ict','droptol',1e-2));
[x4,fl4,rr4,it4,rv4] = pcg(A,b,tol,maxit,L4,L4');
L5 = ichol(A, struct('type','ict','droptol',1e-3));
[x5,fl5,rr5,it5,rv5] = pcg(A,b,tol,maxit,L5,L5'); 

figure
semilogy(0:maxit,rv0./norm(b),'b-','linewidth',2);
hold on
semilogy(0:it3,rv3./norm(b),'b-.','linewidth',2);
semilogy(0:it4,rv4./norm(b),'b--','linewidth',2);
semilogy(0:it5,rv5./norm(b),'b:','linewidth',2);
legend('No Preconditioner','ICT(1e-1)','ICT(1e-2)', ...
   'ICT(1e-3)','Location','SouthEast');

% Note the incomplete Cholesky preconditioner constructed with drop tolerance 1e-2 is denoted as ICT(1e-2).
% As with the zero-fill incomplete Cholesky, the threshold dropping factorization can benefit from modification (i.e. opts.michol = 'on') since the matrix arises from an elliptic partial differential equation. As with MIC(0), the modified threshold based dropping incomplete Cholesky will preserve the action of the preconditioner on constant vectors, that is norm(A*e-L*(L'*e)) will be approximately zero.
% Copyright 2015 The MathWorks, Inc.