%% Using the diagcomp Option
% This example illustrates the use of the diagcomp option of ichol.
% Incomplete Cholesky factorizations of positive definite matrices do not always exist. The following code constructs a random symmetric positive definite matrix and attempts to solve a linear system using pcg.

S = rng('default');
A = sprandsym(1000,1e-2,1e-4,1);
rng(S);
b = full(sum(A,2));
[x0,fl0,rr0,it0,rv0] = pcg(A,b,1e-6,100);

% Since convergence is not attained, try to construct an incomplete Cholesky preconditioner.
% L = ichol(A,struct('type','ict','droptol',1e-3));
% Error using ichol
% Encountered nonpositive pivot.
%
% If ichol breaks down as above, you can use the diagcomp option to construct a shifted incomplete Cholesky factorization. That is, instead of constructing L such that L*L' approximates A, ichol with diagonal compensation constructs L such that L*L' approximates M = A + alpha*diag(diag(A)) without explicitly forming M. As incomplete factorizations always exist for diagonally dominant matrices, alpha can be found to make M diagonally dominant.

alpha = max(sum(abs(A),2)./diag(A))-2

%

L1 = ichol(A, struct('type','ict','droptol',1e-3,'diagcomp',alpha));
[x1,fl1,rr1,it1,rv1] = pcg(A,b,1e-6,100,L1,L1');

%
% Here, pcg still fails to converge to the desired tolerance within the desired number of iterations, but as the plot below shows, convergence is better for pcg with this preconditioner than with no preconditioner. Choosing a smaller alpha may help. With some experimentation, we can settle on an appropriate value for alpha.

alpha = .1;
L2 = ichol(A, struct('type','ict','droptol',1e-3,'diagcomp',alpha));
[x2,fl2,rr2,it2,rv2] = pcg(A,b,1e-6,100,L2,L2');

%
% Now, pcg converges and a plot can show the convergence histories with each preconditioner.

semilogy(0:100,rv0./norm(b),'b.');
hold on;
semilogy(0:100,rv1./norm(b),'r.');
semilogy(0:it2,rv2./norm(b),'k.');
legend('No Preconditioner','\alpha \approx 63','\alpha = .1');
xlabel('Iteration Number');
ylabel('Relative Residual');

% Copyright 2015 The MathWorks, Inc.