function [x, resvec] = cgsolve(A, b, niters)
% cgsolve : Solve Ax=b by conjugate gradients
%
% [x, resvec] = cgsolve(A, b, niters);
%
% Given symmetric positive definite sparse matrix A and vector b, 
% this runs "niters" iterations of 
% conjugate gradient to solve for x in A*x=b.
%
% Optionally, [x, resvec] = ... also returns a vector of 
% the norm of the residual b - A*x at each iteration.
%
%
% The code below follows the slides on the CS111 web site.
%
% John R. Gilbert    revised 13 Jan 2016

n = length(b);
normb = norm(b);

% Initial guess: x = 0
x = zeros(n,1);

% Initial residual: r = b - A*0 = b
r = b;
 
% Initial step is in direction of residual.
d = r;

rtr = r'*r;
resvec = zeros(1,niters);

% Iterate
for iter = 1 : niters
    Ad = A*d;
    alpha = rtr / (d'*Ad);  % Length of step
    x = x + alpha * d;      % Update x to new x
    r = r - alpha * Ad;     % Update r to new residual
    resvec(iter) = norm(r);
    rtrold = rtr;
    rtr = r'*r;
    beta = rtr / rtrold;   
    d = r + beta * d;       % Update d to new step direction
end;
