% PCG Demo
clear all
close all
clc
scrnsz = get(0, 'Screensize');
m = 1000;
for i=1:m
    d(i) = 0.5 + sqrt(i);
end
A = diag(ones(m-100,1),-100) + diag(ones(m-1,1),-1) + diag(d) + ...
    diag(ones(m-1,1),1) + diag(ones(m-100,1),100);
b = ones(m,1);
% without preconditioning
%n = nnz(A);
%kappa = cond(A);
for iter = 1 : 40
    [x, flag, relres] = pcg(A,b,1e-16,iter);
    r(iter) = relres;
end
figure('Position', [5 10 scrnsz(3)-5 scrnsz(4)-60])
semilogy(r,'g+','Linewidth',2)
xlabel('Number of PCG iterations','Fontsize',14);
ylabel('Relative residual','Fontsize',14);
title('Convergence of Diagonally Preconditioned CG Algorithm (accuracy 1e-16)','Fontsize',18);
%t1 = [sprintf('Without preconditioning, nz = %5d, ', n), '\kappa', sprintf(' = %3.2e',kappa)];
t1 = [sprintf('Without preconditioning')];
legend(t1, 3)
drawnow;
% preconditioned
for iter = 1 : 40
    [x, flag, relres] = pcg(A,b,1e-16,iter,diag(diag(A)));
    r(iter) = relres;
end
hold on;
semilogy(r,'ro','Linewidth',2)
%t2 = [sprintf('With preconditioning, nz = %5d, ', n), '\kappa', sprintf(' = %3.2e',kappa)];
t2 = [sprintf('With preconditioning')];
legend(t1, t2, 3)
drawnow;
