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.  In the homework, the matvec will use
% the 5-point finite difference matrix for the 2D model problem.
%
% John R. Gilbert    12 Feb 2006, modified 11 April 2009

n = length(b);
maxiters = max(100,sqrt(n));

% Note: for initial debugging of the routines you might
%       want to set maxiters to something small like 10.

normb = norm(b);
x = zeros(n,1);
r = b;
rtr = r'*r;
d = r;
niters = 0;
while sqrt(rtr)/normb > 1e-4  &&  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;
