function [x,f,g,H] = newmod(func,x0)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Prototype: [x,f,g,H] = newmod(func,x0)
%
% Purpose:   Finds a point satisfying the first-order optimality conditions
%            for a nonlinear function f(x) using modified Newton with Armijo
%            backtracking.  Stopping criteria is ||g||_2 < eps.
%
% 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,H] = func(x)
%
% Output:    x     = optimal point
%            f,g,H = f,g,H at optimal point
%
% Author:  M. Holst (implementation follows P. Gill's notes closely)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
tol     = sqrt(eps);    % tolerance for ||g||_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;              
alpha   = 0;
dmin    = 0;
x       = x0;   
[f,g,H] = feval(func,x);
[V,D]   = eig(H);
D       = diag(D);
dmin    = min(D);         
fprintf(' iter alpha         f(x)      ||g||_2    Min eig\n' )
fprintf(' %3g %9.2e %14.7e %9.2e %9.2e\n', iter, alpha, f, norm(g), dmin);
while ( (norm(g) >= tol) & (iter < itmax) )
   D     = max( abs(D), delta );    
   B     = V*diag(D)*V';
   p     = - B\g;
   gtp   = g'*p;
   alpha = 1;
   xnew  = x + alpha*p;
   [fnew,gnew,Hnew] = feval(func,xnew);
   j = 1;
   while ( (fnew > f + alpha*mu*gtp) & (j <= dpmax) )
      alpha = alpha/2;
      xnew  = x + alpha*p;
      [fnew,gnew,Hnew] = feval(func,xnew);
      j = j + 1;
   end
   x       = xnew;  
   [f,g,H] = feval(func,x);
   [V,D]   = eig(H);
   D       = diag(D);
   dmin    = min(D);         
   iter    = iter + 1;
   fprintf(' %3g %9.2e %14.7e %9.2e %9.2e\n', iter, alpha, f, norm(g), dmin);
end
