  function [x,nrmHist] = rnsteep(A, b, x, nsteps,tol)
% function [x,nrmHist] = rnsteep(A, b, x, nsteps,tol)
% Res. Norm steepest descent 
% function [x,r] = rnsteep(A, x, b, nsteps) 
% A = matrix = 
% x = initial guess on entry -- app sol on return 
% b = rhs, 
% nsteps = number of steps -- r = b - A * x;
 r = b-A*x;
 ro = norm(r);
 nrmHist(1) = ro;
 if (nargin<5), tol = eps;, end
 tol_rel = tol*ro;
%%-------------------- loop
 for i=1:nsteps
	 d = A' * r;
         Ad = A*d; 
	 alp = d'*d / ( Ad' * Ad);
	 x = x + alp * d;
	 r = r - alp* Ad;
         ro = norm(r);
%%-------------------- save norms for plotting
         nrmHist(i+1) = ro;
         if (ro < tol_rel), break, end          
 end 
          
