function [Xpt,h1] = SteepestDescentF(fcn,grad,Xpt,itMax,tol,varargin)
% SteepestDescentF.m
% Function using the speepest descent method for multi-parameter optimization.
% Line search is implemented using the Armijo method.
% fcn       string of a file or file handle of a (scalar) cost function 
% grad      string of a file or file handle of the (vector) gradiant of fcn
% Xpt       parameter vector (initialized to a starting vector)
% itMax     maximum number of iterations allowed
% tol       tolerance goal for the norm of the parameter vector
% varargin  cell array (struct) containing any optional input parameters to
%           be supplied to the cost and the gradiant functions
% OUTPUT:
% Xpt       parameter vector at a minimum
% h1        function at the minimum
% K. Ming Leung, 10/09/04

absErr = 1;     % absolute error in the solution vector
it = 0;         % iteration count
h1 = feval(fcn,Xpt,varargin);

disp('                parameter vector             function value  ');
%fprintf(1,'%10.6f %10.6f %10.6f %10.6f\n',Xpt.',h1);
disp([transpose(Xpt) h1]);  % display parameters and function value

while (absErr > tol && it < itMax)
  it = it +1;
  Step = -feval(grad,Xpt,varargin);    % steepest descent direction
  alpha3 = 1;                   % initial value for alpha3
  h3 = feval(fcn,Xpt+alpha3*Step,varargin);    % function value at alpha3
  while (h3 >= h1)
    alpha3 = 0.5*alpha3;                % keep on halfing alpha3 till h3 < h1  
    h3 = feval(fcn,Xpt+alpha3*Step,varargin);
  end
  alpha2 = 0.5*alpha3;                  % let alpha2 = alpha3/2
  h2 = feval(fcn,Xpt+alpha2*Step,varargin);
  c2 = (h2-h1)/alpha2;                  % coefficient from interpolation
  c22 = (h3-h1)/alpha3;                 % coefficient from interpolation
  c3 = (c22-c2)/(alpha3-alpha2);        % coefficient from interpolation
  alphaTemp = 0.5*(alpha2-c2/c3);       % min of the parabolic interpolant
  
  % The following may be needed when Xpt is far from the min
  if alphaTemp > 0
    alphaMin = alphaTemp;
  else
    alphaMin = alpha3;
  end
  
  DeltaXpt = alphaMin*Step; % change in the parameter vector
  absErr = norm(DeltaXpt);  % norm of the change in the parameter vector
  Xpt = Xpt+alphaMin*Step;  % new parameter vector
  h1 = feval(fcn,Xpt,varargin);    % new function value there
  %fprintf(1,'%10.6f %10.6f %10.6f %10.6f %10.6f\n',Xpt.',h1);
  disp([transpose(Xpt) h1]);  % display parameters and function value
  
end % iteration loop


