function [xrec,frec,drec,srec,hitbound] = gradsampfixed(pars, x0, samprad, ...
    dnormtol, options)
% gradient sampling minimization with fixed sampling radius
%  pars:   problem parameters needed, including pars.m,
%    pars.fgname: name of mfile coding function and gradient
%    pars.fcomparename: name of mfile comparing function value with given value
%  x0: starting point
%  samprad: perturbation parameter for evaluating nearby gradients
%  dnormtol: quitting tolerance for norm(d)
%  options: options defining the algorithm, including:
%    options.maxit: max number of iterations
%    options.xbound: quit if inf-norm(x) hits this
%    options.ngrad: number of gradients per bundle
%    options.prtlevel: print level
%    options.subprob: 'Q' for solving QP subproblem 
%                     'L' for solving LP subproblem 
%
%  Written by M. Overton (overton@cs.nyu.edu), last revised November 2004
fgname = pars.fgname;
prtlevel = options.prtlevel;
ngrad = options.ngrad;
if prtlevel > 0
   fprintf('sampling radius = %7.1e  ',samprad);
end
x = x0;
xrec = [];
frec = [];
drec = [];
srec = [];
hitbound = 0;
maxit = options.maxit;
ngrad = options.ngrad;
subprob = options.subprob;
xbound = options.xbound;
for k=1:maxit
   xrec(:,k) = x;
   % evaluate function and gradients (including function and gradient at x)
   [fbundle, gbundle] = getbundle(x, samprad, ngrad, pars);   
   f = fbundle(1);     % function values at perturbed points are discarded
   frec(k) = f;
   if subprob == 'Q'   % solve QP subproblem
      [w,d] = qpsubprob(gbundle, prtlevel); % recommended
   else                % solve LP subproblem
      [w,d] = lpsubprob(gbundle, prtlevel); % not recommended
   end
   g = gbundle(:,1)'*d;  % gradient value at current point, without perturbation
   dnorm = norm(d);      % will normalize d below
   drec(k) = dnorm;
   if dnorm < dnormtol    
      if prtlevel > 0
         fprintf('success in %d iterations: norm(d) = %5.1e \n', ...
         k, norm(d));
      end
      step = 0; 
   elseif g >= 0 
      if prtlevel > 0
         fprintf('quit at iteration %d: d not descent direction: g*d = %e \n', ...
         k, g);
      end
      step = 0; 
   elseif isnan(g)
      fprintf('quit at iteration %d: g is NaN', k)
      step = 0; 
   else 
      % nondifferentiable line search - normalize d first
      % important: otherwise huge steps can make trouble
      % (makes a big difference for probdef_exp(3), for example)
      % small steps can also cause difficulties, especially if the line
      % search does not have a doubling feature, so we always normalize d
      d = d/dnorm; 
      smax = stepmax(x, d, xbound);
%
% This line search does not compute function values.
% All it does is find a point x + step*d for which the new function value
% is lower than the current function value f.  In many applications, it
% is much cheaper to generate bounds on the function value than to compute
% the function value exactly, and since the line search may have to try
% MANY steplengths till it finds an acceptable one, a lot of time can be
% saved.  The actual function value at the new point is computed along with
% the gradient bundle at the top of the loop.
%
      [step, fail] = ndlinesch(x, f, d, smax, pars);  % formerly ndlinesch2
      x = x + step*d;
      if fail == 1   % fail = -1 means didn't bracket min but that's ok
                     % fail = -2 means stepmax was reached
         if prtlevel > 0
            fprintf('line search did not get reduction at iteration %d \n', k);
         end
      end
   end
   srec(k) = step;
   if step == 0 % quit if d small or not descent direction or line search failed
      return
   elseif step == smax % quit if hit bound on x
      disp('hit bound on x')
      hitbound = 1;
      xrec(:,k+1) = x;
      frec(k+1) = feval(fgname, x, pars); % computed bound in line sch
      return    
   end
end
xrec(:,maxit + 1) = x;    % final iterate.  Starting point is in 1st position.
frec(maxit + 1) = feval(fgname, x, pars); % computed bound in line sch
if prtlevel > 0
   fprintf('max iterations exceeded, norm(d) = %g \n', dnorm);
end
