function [xn,ro,fxn,info]=goldstein(x,d,fx,gfx,r,f)
%
% function [xn,ro,fxn,info]=goldstein(x,d,fx,gfx,r,f)
%
%     Copyright: Eduardo Casas (Octubre, 2000)
%
% Esta funcion implementa la regla de Goldstein. Utiliza
% como funcion auxiliar: nuevot.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.
% r      Reduccion prevista.
% f      Funcion a minimizar.
%
% 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.
% info   Variable que toma el valor 0 si [xn,ro,fxn]
%        ha sido calculado en el numero maximo de iteraciones
%        previsto y -1 en caso contrario.
%
maxiter=20;
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;
ro=min([max([a,-2*r/gfxd]),b]);
for i=1:maxiter
   xn=x+ro*d;
   fxn=feval(f,xn);
   if fxn > fx + sigma*ro*gfxd
      b=ro;
      info=2;
   elseif fxn < fx + beta*ro*gfxd;
      a=ro;
   else
      info=0;
      return;
   end
   if i == 1
      h=[fx;gfxd;fxn];
      t=ro;
      ro=nuevot(a,b,t,h,info);
      t=[0;0;t];
      h(2)=fx;
   else
      t=[t(2);t(3);ro];
      h=[h(2);h(3);fxn];
      ro=nuevot(a,b,t,h,info);
   end
end
info=-1;
%
%
function tn=nuevot(a,b,t,h,info)
%
% function tn=nuevot(a,b,t,h,info)
%
% Esta funcion calcula el minimo relativo tn en (a,b) de un
% polinomio de grado 2. Si t es un escalar, entonces se supone
% que h(1) representa el valor del polinomio en 0, h(2) su
% derivada en el mismo punto y h(3) el valor del polinomio en
% t. Si t es un vector de tres componentes, entonces h es un
% vector que contiene los valores del polinomio en los nodos de
% t:         p(t(i)) = h(i).
% Si tal minimo no exise se tomara el valor (a+b)/2 o el minimo
% de b y 10a dependiendo de si info=2 o 1.
%
n=length(t);
if n==1
   c=[(h(3)-h(2)*t-h(1))/t^2;h(2);h(1)];
else
   c=[t.^2,t,ones(3,1)]\h;
end
if c(1) > eps
   tn=-c(2)/(2*c(1));
   if info == 1
      tn=min([max([tn;10*a]);b]);
   else
      r=(b-a)/10;
      tn=min([b-r;max([a+r,tn])]);
   end
else
   if info == 1
      tn = min([b,10*a]);
   else
      tn = a+(b-a)/2;
   end
end
