function [xnew,gradnew,alpha,p] = bfgs(searchmtd,fnhess,fngrad,fn,gradold,xold)
% usage: [xnew,gradnew,alpha,p] = bfgs(searchmtd,fnhess,fngrad,fn,gradold,xold)
% description: does a single step of BFGS method with input parameters:
% searchmtd: line search method
% fnhess: function hessian name
% fngrad: function grad name
% fn: function name
% gradold: gradient at starting point
% xold: starting point

% global variables:
global parms;
global H; % inverse of current Hessian approx

% local variables:
% xold,xnew: old,new x values
% gradold,gradnew: gradf(x*)
% alpha: step length used
% p: descent direction used

p = -H*gradold;
n = length(xold);
if (p'*gradold >= 0) % not a descent direction?
  disp('Restart due to non-descent direction')
  H = eye(n);
  p = -gradold;
end 
alpha = feval(searchmtd,fngrad,fn,xold,p);
xnew = xold + alpha*p;
gradnew = feval(fngrad,xnew);
s = xnew-xold;
y = gradnew - gradold;
rho = y'*s;
% update H for next iteration
if (rho <= 0) % rho negative? shouldn't happen
  disp('Did a matrix restart due to negative rho')
  H = eye(n); % do a restart
else
  mtmp = eye(n)-s*y'/rho;
  H = mtmp*H*mtmp + s*s'/rho;
end