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

%  function [u,rho] = steepest_descents(A,u0,f,relerrtol,maxiter)
%  Calculate solution to Au = f via steepest descents, 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 convergenc
      return
    end
    w = A*r;
    s = rho(iter)/(r'*w);
    u = u + s*r;
    r = r - s*w;          % Shortcut for calculating residual 
  end

%  Error message

  ['Steepest descents did not converge after ' num2str(maxiter) ' iterations']
