[x,f,g] = steepest(func,x0)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Prototype: [x,f,g] = steepest(func,x0)
%
% Purpose:   Performs steepest descent to find a minimizer of a given function.
%
% Input:     x0   = initial guess at the minimizer.
%            func = string containing the name of a m-file that
%                   defines the function f, with prototype: [f,g] = func(x)
%
% Output:    x   = optimal point
%            f,g = f,g at optimal point
%
% Author:  M. Holst (implementation follows P. Gill's notes closely)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
tol   = sqrt(eps);    % tolerance for ||F||_2
itmax = 50;           % limit on function evaluations
dpmax = 20;           % limit on damping steps
mu    = 1/4;          % sufficient decrease paramter
delta = 1.0e-10;      % smallest eigenvalue cut-off
iter  = 0;              
nf    = 0;
alpha = 0;
x     = x0;   
[f,g] = feval(func,x);
clc                                 
fprintf( ' Itn   Nf      Step      f(x)         Norm g\n' )
str1 = sprintf ( ' %3g %4g %9.2e', iter, nf, alpha );
str2 = sprintf ( ' %14.7e %9.2e' , f, norm( g ) );
disp( [str1 str2 ] )                
while ( (norm(g) >= tol) & (iter < itmax) )
   p           = - g;
   gtp         = g'*p;
   alpha       = 1;
   xnew        = x + alpha*p;
   [fnew,gnew] = feval(func,xnew);
   jf          = 1;
   while  ( (fnew > f+alpha*mu*gtp) & (jf < dpmax) )
   % while  ( (abs(fnew) > (1-mu*alpha)*abs(f)) & (jf < dpmax) )
      alpha       = alpha/2;
      xnew        = x  + alpha*p;
      [fnew,gnew] = feval(func,xnew);
      jf          = jf + 1;
   end
   nf    = nf + jf;
   x     = xnew;  
   [f,g] = feval(func,x);
   iter  = iter + 1;
   str1  = sprintf ( ' %3g %4g %9.2e', iter, nf, alpha );
   str2  = sprintf ( ' %14.7e %9.2e', f, norm( g ) );
   disp( [str1 str2 ] )
end
if ( norm(g) <= tol )
   fprintf( '\nOptimal at iteration %g.  ', iter )
   fprintf( '  Objective = %14.7e\n\n', f );
else
   fprintf( '\nTerminated.  Iteration or backtrack limit exceeded  \n\n' )
end

