%% Example of linear CG method for a matrix with clustered eigenvalues

clc, clear all, clf

% This is the implementation of the 
% conjugate gradient method (Algorithm 5.1 in
% Nocedal and Wright "Numerical Optimization") 
% for solving linear systems Ax=b, with A symmetric and 
% positive definite matrix, 

Centers=[1 10 30];
Nclusters=length(Centers);

A=[];
for i=1:Nclusters
    B=rand(10); B=B'*B;
   
    B=(B/max(eig(B))+10*eye(10))/10+Centers(i)*eye(10);
     eig(B)
    A=blkdiag(A,B);
end

b=rand(Nclusters*10,1);

n=length(b);
xstar = A\b; % uses the 'exact' solution for computing the error 
% done at each iteration

% choose the number of iterations to be performed
Niter = n;

res = zeros(Niter,1);
err = zeros(Niter,1);
errA = zeros(Niter,1);


x = zeros(n,1); % initial guess
r = A*x-b;

p = -r; % initial direction (steepest descent)
for k = 1:Niter
    alpha = -(r'*p)/(p'*A*p);
    x = x + alpha*p;
    r = A*x-b;
    res(k) = r'*A*r; % computes the A-norm of the residuals
   % err(k) = norm(xstar-x); % computes the errors
    errA(k) = (x-xstar)'*r; % computes the error in the A-norm
    beta = (r'*A*p)/(p'*A*p);
    p = -r+beta*p;  
end
xCG = x;

subplot(2,1,1)
plot(eig(A),zeros(n,1),'o')
grid on

% subplot(3,1,2)
% semilogy(res,'-'); grid on
% axis([0 n 10^(-16) max(res)*1.5])
% title('Log Plot of the residual:  $\| r^k\|_A^2$ ', 'interpreter','latex') 


subplot(2,1,2)
semilogy(errA,'-o'); grid on
xlabel('iteration k')
axis([0 n 10^(-16) max(errA)*1.5])
title('Log Plot of the error:  $\| x^k-x^*\|_A^2$','interpreter','latex') 
