function cg
% CG on Ax = b

% Define matrix A
dim = 1000; % matrix dimension
%A = hilbertMatrix(dim);
A = generateMatrix(dim, 1);

b = ones(dim, 1);

ep = 1.0E-6;
x0 = zeros(dim,1);
x = cg_loop(A,b,x0,ep);

end


function x = cg_loop(A,b,x,ep)
x1 = A\b; % "Exact" solution
fprintf('#it\t||r||\t\tlog(||x-x*||_A)\n')

r = A*x - b;
p = -r;
rr= r'*r;
k = 0;
while norm(r) > ep
    Ap    = A*p;
    alpha = rr/(p'*Ap);
    x     = x+alpha*p;
    r     = r+alpha*Ap;
    rr1   = r'*r;
    beta  = rr1/rr;  
    p     = -r+beta*p;
    rr    = rr1;
    k     = k+1;
    lognorm(k) = log(sqrt((x1-x)'*A*(x1-x)));
    fprintf('%d:\t%8f\t%8f\n', k, norm(r), lognorm(k));
end

fprintf('CG converged after %d iterations\n', k);
fprintf('Relative error compared to direct solve: %e\n', norm(x1-x)/norm(x1));
plot(1:k,lognorm)

end


function A = hilbertMatrix(dim)
% A_ij = 1/(i+j-1)
A = zeros(dim);
for i=1:dim
    for j=1:dim
        A(i,j) = 1/(i+j-1);
    end
end
fprintf('Condition number of A: %g\n', cond(A));

end


function A = generateMatrix(n, type)
% Create matrix of dimension n (multiple of 5) and cond(A)=100
% Matrix types:
%   1: 5 clustered eigenvalues
%   2: 5 clustered and small perturbation
%   3: uniformly distributed eigenvalues
if (mod(n,5) ~= 0) error('Dimension should be multiple of 5'); end 
rng('default'); % To produce same matrix every time
M = rand(n); % Random nxn matrix with elements in the range [0,10]
[Q, R] = qr(M); % Q is now unitary

if (type == 1) || (type == 2)
    rep = n/5;
    e = ones(rep,1);
    D = [1*e; 2*e; 10*e; 20*e; 100*e];
    if (type == 2)
        for i=1:n
            D(i) = D(i) + 0.001*(0.5-rand())*D(i);
        end
    end
elseif (type == 3)
    D = 1:(99/(n-1)):100;
else 
    error('Type should be 1,2 or 3')
end      
D = diag(D);
A = Q'*D*Q;

end