% Vicol Vlad, Project 5, Problem 5; small changes by M.O.
% gradient method

% this function solves a linear system Ax=b using the gradient method.
% The parameters are A, b, and a starting vector.
% The function returns the iteration history (the 2-norm of the residual
% at each step), the number of iterations needed, and the solution x.

function [resid, k, xk] = gradientmeth (A,b,st)

  n = length(b);

  xk = st;
  rk = b;
  ak = 0;

  iterlimit = 50; % Maximum number of iterations
  err=10^(-9);     % Desired tolerance
  resid=[];

  k=1;
  while ( (norm(rk) >= err) & (k<=iterlimit) )   
    % while the residual is not small enough and we are not diverging
    rk = b - A*xk;
    resid(k) = norm(rk);
    ak = (rk'*rk) / (rk'*A*rk);
    xk = xk + ak*rk; 
          
    k = k+1;
  end
  k = k-1;

end

