%  Fletcher-Reeves (hopeless!) and Polak-Ribierre (less so)

clear all

x = -1.2 :.01 : 1.2;
y = -0.5 :.01 : 1.5;
[X,Y]=meshgrid(x,y);
Z = l24f(X,Y);
cvals = -1:0.5:10;
contour(x,y,Z,cvals)

hold on

tau = 0.5;
beta = 0.1;
xk = -1.2;
yk = 1.0;
xi(1) = xk;
yi(1) = yk;

 restart = 0 ;
 for k = 2:100;
  restart = restart + 1 ;
  fk = l24f(xk,yk);
  gxk = -40*yk*xk + 40*xk^3 -2 + 2*xk;
  gyk = 20*yk - 20*xk^2 ;
  g = [gxk,gyk];
  if restart > 20
    restart = 1 ;
  end
  if restart == 1
    gtg = gxk^2 + gyk^2 ;
    p = - g ;
  else
    oldgtg = gtg ;
    gtg = gxk^2 + gyk^2 ;
    gtgm1 = gxk * gxkold + gyk * gykold ;
%   beta = gtg / oldgtg ;                          % Fletcher-Reeves
    beta = ( gtg - gtgm1 ) / oldgtg ;              % Polak-Ribierre
    p = - g + beta * p ;
  end  
  pxk = p(1);
  pyk = p(2);
  gtp = pxk*gxk + pyk*gyk ;
  if gtp > 0
    pxk = - p(1);
    pyk = - p(2);
    gtp = pxk*gxk + pyk*gyk ;
  end
  alpha = 1;
  while l24f(xk+alpha*pxk,yk+alpha*pyk) > fk + beta*alpha*gtp
    alpha = alpha * tau;
  end
  alpha ;
  gxkold = gxk ;
  gykold = gyk ;
  xk = xk+alpha*pxk ;
  yk = yk+alpha*pyk ;
  xi(k) = xk;
  yi(k) = yk;
 end


gtg

plot(xi,yi,'-r.')
xsol=1;
ysol=1;
plot(xsol,ysol,'ko')

hold off

print -depsc eps/l25
