function x = gradit(A,b,x,rtol,atol,maxit)
r = b-A*x; % residual = search direction
for k=1:maxit
  p = A*r;
  ts = (r'*r)/(r'*p); % optimal line search parameter
  x = x + ts*r;       % update approximate solution 
  cn = (abs(ts)*norm(r); % norm of correction
  if ((cn < rtol*norm(x)) || (cn < atol)), return; end
  r = r - ts*p; % update residual
end
