function [u,rho] = CG(A,u0,f,relerrtol,maxiter)

%  function [u,rho] = CG(A,u0,f,relerrtol,maxiter)
%  Calculate solution to Au = f via conjugate gradients, where A is spd, from
%  an initial guess u0.
%  Other inputs:
%    relerrtol   u assumed convergent once norm(f - Au)/norm(f) < relerrtol
%    maxiter     Maximum number of iterations
%  Other outputs:
%    rho  Sequence of squared residuals (length maxiter, 0 for indices 
%          greater than the number of iterations required for convergence.

  u = u0;
  tol = (relerrtol^2)*(f'*f);
  rho = zeros(maxiter,1);
  r = f - A*u;

  for iter = 1:maxiter   
    rho(iter) = r'*r;   % rho is squared residual before iteration number iter.
    if (rho(iter) < tol)  % Test for convergence
      return
    end
    if iter == 1               % On first step, search dir is down gradient.
      v = r;
    else                       
      beta = rho(iter)/rho(iter-1);
      v = r + beta*v;          % Update search direction
    end
    w = A*v;
    s = rho(iter)/(v'*w);
    u = u + s*v; 
    r = r - s*w;          % Shortcut for calculating residual
  end

%  Error message

  ['Conjugate gradients did not converge after ' num2str(maxiter) ' iterations']
