function [xn,ro,fxn,gfxn,info]=wolfe(x,d,fx,gfx,ro,f,gf)
%
% function [xn,ro,fxn,gfxn,info]=wolfe(x,d,fx,gfx,ro,f,gf)
%
%     Copyright: Eduardo Casas (Octubre, 2000)
%
% Esta funcion implementa la regla de Wolfe. Utiliza
% como funcion auxiliar: intcub.m. El significado de las
% variables es el siguiente:
%
% Variables de Entrada:
%
% x      Punto actual de las iteraciones.
% d      Direccion de descenso.
% fx     Valor de la funcion f en el punto x.
% gfx    Gradiente de la funcion f en el punto x.
% ro     Punto en el que se inician las iteraciones.
% f      Funcion a minimizar.
% gf     Funcion que calcula el gradiente de f.
%
% Variables de Salida:
%
% xn     Nuevo punto calculado: xn = x + ro*d.
% ro     Valor del paso calculado.
% fxn    Valor de la funcion f en el nuevo punto xn.
% gfxn   Valor de la funcion gf en el nuevo punto xn.
% info   Variable que toma el valor 0 si [xn,ro,fxn,gfxn]
%        ha sido calculado en el numero maximo de iteraciones
%        previsto y -1 en caso contrario.
%
maxiter=50;
sigma=0.001;
beta=0.8;
info=1;
nd=norm(d);
a=max([eps*norm(x)/nd,eps]);
b=sqrt(realmax)/max([nd,1]);
gfxd=gfx'*d;
s1=sigma*gfxd;
s2=beta*gfxd;
ro=min([max([a,ro]),b]);
t=[0;0];
h=[fx;gfxd;fx;gfxd];
for niter=1:maxiter
   xn=x+ro*d;
   fxn=feval(f,xn);
   gfxn=feval(gf,xn);
   gfxnd=gfxn'*d;
   if fxn > fx + ro*s1
      b=ro;
      info=2;
   elseif gfxnd < s2;
      a=ro;
   else
      info=0;
      return;
   end
   t=[t(2),ro];
   h=[h(3);h(4);fxn;gfxnd];
   ro=intcub(t,h,a,b,info);
end
info=-1;
%
%
function ro=intcub(t,h,a,b,info)
%
% function ro=intcub(t,h,a,b,info)
%
%
% Esta funcion calcula el minimo relativo ro en (a,b) de un
% polinomio de grado 3 que en los nodos t(1) y t(2) toma los
% valores h(1) y h(3) y su derivada en dichos puntos es h(2)
% y h(4) respectivamente. Si tal minimo no existe se tomara el
% valor (a+b)/2 o el minimo de b y 10a dependiendo de si info=2
% o 1.
%
l=t(2)-t(1);
x=2*(h(1)-h(3))+(h(2)+h(4))*l;
y=3*(h(3)-h(1))-(2*h(2)+h(4))*l;
z=h(2)*l;
d=y^2-3*x*z;
if d <= 0
   ind=1;
else
   r=abs(y)+sqrt(d);
   if y > 0
      ro=t(1)-l*z/r;
      ind=0;
   elseif y<= 0 & abs(x) >= eps
      ro=t(1)+l*r/(3*x);
      ind=0;
   else
      ind=1;
   end
end
if ind == 0
   if info == 1
      ro=min([max([ro,10*a]),b]);
   else
      r=(b-a)/10;
      ro=min([max([ro,a+r]),b-r]);
   end
else
   if info == 1
      ro=min([10*a,b]);
   else
      ro=a+(b-a)/2;
   end
end
