% CG Demo
clear all
close all
clc
scrnsz = get(0, 'Screensize');
m = 500;
A = triu(-1+2*randn(m,m),1);
A = A + A';
A = A + eye(m);
b = rand(m,1);
% tau = 0.2
tau = 0.2;
A(find(abs(A)>tau))=0;
A = A + eye(m);
n = nnz(A);
kappa = cond(A);
for iter = 1 : 20
    [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 CG iterations','Fontsize',14);
ylabel('Relative residual','Fontsize',14);
title('Convergence of Conjugate Gradient Algorithm (accuracy 1e-16)','Fontsize',18);
t1 = ['\tau',sprintf(' = %3.2f, nz = %5d, ', tau,n), '\kappa', sprintf(' = %3.2e',kappa)];
legend(t1, 3)
drawnow;
% tau = 0.1
tau = 0.1;
A(find(abs(A)>tau))=0;
A = A + eye(m);
n = nnz(A);
kappa = cond(A);
for iter = 1 : 20
    [x, flag, relres] = pcg(A,b,1e-16,iter);
    r(iter) = relres;
end
hold on;
semilogy(r,'ro','Linewidth',2)
t2 = ['\tau',sprintf(' = %3.2f, nz = %5d, ', tau,n), '\kappa', sprintf(' = %3.2e',kappa)];
legend(t1, t2, 3)
drawnow;
% tau = 0.05
tau = 0.05;
A(find(abs(A)>tau))=0;
A = A + eye(m);
n = nnz(A);
kappa = cond(A);
for iter = 1 : 20
    [x, flag, relres] = pcg(A,b,1e-16,iter);
    r(iter) = relres;
end
hold on;
semilogy(r,'bx','Linewidth',2)
t3 = ['\tau',sprintf(' = %3.2f, nz = %5d, ', tau,n), '\kappa', sprintf(' = %3.2e',kappa)];
legend(t1, t2, t3, 3)
drawnow;
% tau = 0.01
tau = 0.01;
A(find(abs(A)>tau))=0;
A = A + eye(m);
n = nnz(A);
kappa = cond(A);
for iter = 1 : 20
    [x, flag, relres] = pcg(A,b,1e-16,iter);
    r(iter) = relres;
end
hold on;
semilogy(r,'m*','Linewidth',2)
t4 = ['\tau',sprintf(' = %3.2f, nz = %5d, ', tau,n), '\kappa', sprintf(' = %3.2e',kappa)];
legend(t1, t2, t3, t4, 3)
drawnow;
% tau = 0.01, preconditioned
for iter = 1 : 20
    [x, flag, relres] = pcg(A,b,1e-16,iter,diag(diag(A)));
    r(iter) = relres;
end
hold on;
semilogy(r,'cd','Linewidth',2)
t5 = ['\tau',sprintf(' = %3.2f, nz = %5d, ', tau,n), '\kappa', sprintf(' = %3.2e',kappa)];
legend(t1, t2, t3, t4, t5, 3)
drawnow;
hold off;