function [xnew,gradnew,alpha,p] = gaussnewton(searchmtd,fnjac,fnvec,fngrad,fn,gradold,xold)
% usage: [xnew,gradnew,alpha,p] = newton(searchmtd,fnjac,fnvec,fngrad,fn,gradold,xold)
% description: does a single step of Gauss-Newton method. Returns
% new x, its gradient, step length and search direction.
% Input parameters:
% searchmtd: line search method
% fnjac: function jacobian name
% fnvec: vector function r(x)
% fnhess: function hessian name
% fngrad: function grad name
% fn: function name for f(x)=1/2||r(x)||^2
% gradold: gradient at starting point
% xold: starting point

% global variables:
global parms;

% local variables:
% xnew: new x value
% gradnew: gradf(xnew)
% alpha: step length used
% p: descent direction used

% this is expensive:
%gnmin = parms.gnmin;
%p = -pinv(feval(fnjac,xold),gnmin)*feval(fnvec,xold);
% cheaper alternative:
p = -feval(fnjac,xold)\feval(fnvec,xold);
if (p'*gradold >= 0) % not a descent direction? so restart
  disp('Restart due to non-descent direction');
  p = -gradold;
end
alpha = feval(searchmtd,fngrad,fn,xold,p);
xnew = xold + alpha*p;
gradnew = feval(fngrad,xnew);