function xnew = dbrent_linesearch (x, d, f, df)

  % Parameters are
  % x:    Current location
  % d:    A descent direction
  % f:    The function to be minimized
  % df:   The gradient of f (as a column vector)
  % xnew: A minimum of f along the direction d

  % We first bracket the minimum

  a = 0;
  b = 1;

  fa = feval (f, x);
  fb = feval (f, x+d);

  if fa > fb
     c  = 2;
     fc = feval (f, x+c*d);
     while fb > fc
       c = 2*c;
       fc = feval (f, x+c*d);
     end
  else
     c  = b;
     fc = fb;
     b  = 0.5;
     fb = feval (f, x+b*d);
     while fb > fa
       b = 0.5*b;
       fb = feval (f, x+b*d);
     end
  endif

  % Now Dekker-Brent with derivative information
  % from Numerical Recipes, see Chapter 10.3

  zeps  = 1e-10;
  tol   = 1e-8;
  maxit = 100;
  alpha = -1;
  prevs = 0;

  v   = w   = b;
  fv  = fw  = fb;
  dfv = dfw = dfb = feval (df, x+b*d)' * d;

  for i = 1:maxit

    m    = (a+c)/2;      % Midpoint of the bracketing interval
    tol1 = tol*b + zeps;
    tol2 = 2*tol1;

    if abs(b-m) <= tol2 - (b-a)/2   % Interval has shrunk enough
       alpha = b;                   % Get out with minimum at b
       break;    
    endif

    if abs(prevs) > tol1     % We have noticably moved the previous step
       s1 = 2*(b-a);
       s2 = s1;

       if dfw != dfb
	  s1 = (w-b)*dfb/(dfb-dfw);  % Secant method between b and w
       endif
       if dfv != dfb
	  s1 = (v-b)*dfb/(dfb-dfv);  % Secant method between b and v
       endif

       u1 = b+s1;  % New point from first secant method
       u2 = b+s2;  % New point from second secant method

       % Accept steps if the remain in bracket, and have right direction
       ok1 = (a-u1)*(u1-c) > 0 && dfb*s1 < 0;
       ok2 = (a-u2)*(u2-c) > 0 && dfb*s2 < 0;

       pprevs = prevs;  % Movement on step before last
       prevs  = s;      % Update current movement as previous movement

       if ok1 || ok2

	  if ok1 && ok2           % When both possible steps are OK
	     if abs(s1) < abs(s2) % Take the smallest of the two
	        s = s1;
             else
	        s = s2;
	     endif
	  elseif ok1              % Take the good step
	     s = s1;
	  else
	     s = s2;
	  endif

	  if abs(s) <= abs(prevs/2)
	     u = b+s;
	     if u-a < tol2 || c-u < tol2
	        s = tol1 * sign(m-b);
             endif
	  endif

       else                       % Resort to bisection
          
	  if dfb >= 0
	     prevs = a-b;
	  else
	     prevs = c-b;
	  endif
	  s = prevs/2;

       endif 

    else % Resort to bisection (same as above)

       if dfb >= 0
	  prevs = a-b;
       else
          prevs = c-b;
       endif
       s = prevs/2;

    endif 

    if abs(s) >= tol1 % We are about to take a step of reasonable length
       u  = b + s;
       fu = feval (f, x + u*d);
    else
       u  = b + tol1*sign(s);
       fu = feval (f, x + u*d);
       if fu > fb    % Miniumum step in downhill direction is uphill
	  alpha = b;
	  break;
       endif
    endif

    dfu = feval (df, x + u*d)' * d;

    if fu <= fb
       if u >= b
	  a = b;
       else
	  c = b;
       endif
       v=w; fv=fw; dfv=dfw;
       w=b; fw=fb; dfw=dfb;
       b=u; fb=fu; dfb=dfu;
    else
       if u < b
	  a = u;
       else
	  c = u;
       endif
       if fu <= fw || w == b
	  v=w; fv=fw; dfv=dfw;
	  w=u; fw=fu; dfw=dfu;
       elseif fu < fv || v == b || v == w
	  v=u; fv=fu; dfv=dfu;
       endif
    endif 

  end

  if alpha >= 0
     xnew = x + alpha*d;
  else
     disp ('No convergence in Dekker-Brent');
  endif

end 