function [x, niters] = cgsolve(A, b)
% cgsolve : Solve Ax=b by conjugate gradients
%
% [x, niters] = cgsolve(A, b);
%
% Given symmetric positive definite sparse matrix A and vector b, 
% this runs conjugate gradient to solve for x in A*x=b.
% It iterates until the residual norm is reduced by 10^-6,
% or for at most max(100,sqrt(n)) iterations.
%
% This is a sequential Matlab template -- CS240A homework 3 is
% to implement this in parallel, including the sparse matvec.
%
% The code below follows the slides from the 8 Feb 2006 class.
%
% John R. Gilbert    12 Feb 2006

n = length(b);
maxiters = max(100,sqrt(n));
normb = norm(b);
x = zeros(n,1);
r = b;
rtr = r'*r;
d = r;
niters = 0;
while sqrt(rtr)/normb > 1e-6  &&  niters < maxiters
    niters = niters+1;
    Ad = A*d;
    alpha = rtr / (d'*Ad);
    x = x + alpha * d;
    r = r - alpha * Ad;
    rtrold = rtr;
    rtr = r'*r;
    beta = rtr / rtrold;
    d = r + beta * d;
end;