function [xnew,gradnew,alpha,p] = broyden(searchmtd,fnjac,fnvec,fngrad,fn,gradold,xold)
% usage: [xnew,gradnew,alpha,p] = broyden(searchmtd,fnjac,fnvec,fngrad,fn,gradold,xold)
% description: does a single step of Broyden 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; % current jacobian approx

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

rxold = feval(fnvec,xold);
p = -H\rxold; % this isn't the most efficient approach
if (p'*gradold >= 0) % not a descent direction?
  disp('Direction change due to non-descent direction')
  p = -p;
end 
alpha = feval(searchmtd,fngrad,fn,xold,p);
xnew = xold + alpha*p;
gradnew = feval(fngrad,xnew);
s = xnew - xold;
y = feval(fnvec,xnew) - rxold;
rho = s'*s;
% update H for next iteration
if (rho <= 0) % rho negative? shouldn't happen
  disp('Did a matrix restart due to negative rho')
  H = feval(fnjac,xold); % do a restart
else
  H = H + (y-H*s)*s'/rho;
end
