%% Incomplete Cholesky Factorization
% This example generates an incomplete Cholesky factorization.
% Start with a symmetric positive definite matrix, A:

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

%
% A is the two-dimensional, five-point discrete negative Laplacian on a 100-by-100 square grid with Dirichlet boundary conditions. The size of A is 98*98 = 9604 (not 10000 as the borders of the grid are used to impose the Dirichlet conditions).
% The no-fill incomplete Cholesky factorization is a factorization which contains only nonzeros in the same position as A contains nonzeros. This factorization is extremely cheap to compute. Although the product L*L' is typically very different from A, the product L*L' will match A on its pattern up to round-off.

L = ichol(A);
norm(A-L*L','fro')./norm(A,'fro')

%

norm(A-(L*L').*spones(A),'fro')./norm(A,'fro')

%
% ichol may also be used to generate incomplete Cholesky factorizations with threshold dropping. As the drop tolerance decreases, the factor tends to get more dense and the product L*L' tends to be a better approximation of A. The following plots show the relative error of the incomplete factorization plotted against the drop tolerance as well as the ratio of the density of the incomplete factors to the density of the complete Cholesky factor.

n = size(A,1);
ntols = 20;
droptol = logspace(-8,0,ntols);
nrm = zeros(1,ntols);
nz = zeros(1,ntols);
nzComplete = nnz(chol(A,'lower'));
for k = 1:ntols
    L = ichol(A,struct('type','ict','droptol',droptol(k)));
    nz(k) = nnz(L);
    nrm(k) = norm(A-L*L','fro')./norm(A,'fro');
end
figure
loglog(droptol,nrm,'LineWidth',2)
title('Drop tolerance vs norm(A-L*L'',''fro'')./norm(A,''fro'')')
figure
semilogx(droptol,nz./nzComplete,'LineWidth',2)
title('Drop tolerance vs fill ratio ichol/chol')

% The relative error is typically on the same order as the drop tolerance, although this is not guaranteed.
% Copyright 2015 The MathWorks, Inc.