function x = cg(evalA,b,x,tol,maxit)
% x supplies initial guess, maxit maximal number of CG steps
% evalA must pass a handle to a MATLAB function realizing A*x
r = b - evalA(x); rho = 1; n0 = norm(r);
for i = 1 : maxit
   rho1 = rho; rho = r' * r; 
   if (i == 1), 
     p = r; 
   else
     beta = rho/rho1; 
     p = r + beta * p;   % update of A-orthogonal basis vectors
   end
   q = evalA(p); 
   alpha = rho/(p' * q);
   x = x + alpha * p;    % update of approximate solution
   if (norm(b-A*x) <= tol*n0) % termination criterion
     return; 
   end 
   r = r - alpha * q;    % update of residual
end  