  function P4A()
% Illustrates pcg with incomplete Cholesky preconditioner obtained via LDLinc.
clc

% Generate a symmetric positive definite Ax = b problem...
GridSize = 60;
Region = {'S','L','C','D','A','H', 'B'}; 
k=1;
A = delsq(numgrid(Region{k},GridSize));
[n,n] = size(A);
b = randn(n,1);
fprintf('Solve Ax = b using pcg with various preconditioners.\n\n')
fprintf('A = delsqr(numgrid(''%1s'',%1d)\n\n',Region{k},GridSize)
fprintf('n = %1d, nnz(A) = %1d\n\n',n,nnz(A))

% Key parameters for the preconditioned congugate gradient method...
% Maximum number of iterations...
itMax = 1000;
% Stop if || Ax - b||/||b| <= tol
tol = .00000000001;

% Without any preconditioner...
G = speye(n);
tic
[x,flag,relres,itLast] = pcg(A,b,tol,itMax,G,G');
t = toc;
fprintf('  Preconditioner    Storage      relres        itLast     \n')
fprintf('------------------------------------------------------------\n')
fprintf(' %15s    %5d     %10.3e        %4d       \n','None',0,relres,itLast)

for p = [0 5 10 20]
   % Incomplete LDL...
   [L,d]  = LDLinc(A,p);
   % Form Cholesky factor because that is what pcg wants...
   L = L*diag(sqrt(d));
   [x,flag,relres,itLast] = pcg(A,b,tol,itMax,L,L');
   s = sprintf('LDLinc(A,%1d)',p);
   fprintf(' %15s    %5d     %10.3e        %4d       \n',s,nnz(L),relres,itLast)  
end


