function [success, iter,fcnt,gcnt,fxnew,ginf ] = BB(prob,x)

% *--------------------------------------------------------------
% * March 14,  2005
% * ===============       
% *
% * Program associated to the paper:
% *
% * A New Gradient Descent Method with an Anticipative Scalar 
% * Approximation of Hessian for Unconstrained Optimization.
% *
% * 
% * This is a draft code !   
% * There is no warranty !
% *--------------------------------------------------------------           
      eps  = 10.d0^(-6);  
      epsf = 10.d0^(-20);    
      
%*-------------------------------------- Parameters in backtracking      
      gama=0.0001d0;        
      s=0.8d0 ; 
	     etak=0.7d0;
	  
	     lmdamax =10.d0^30;
	     lmdamin =10.d0^(-30);	  
            
      delta=10.d0  ;
      maxiter=10000 ; 
      maxfgcnt=20000	;  

%*---------------------------------------------  Here Start Experiments                        

        itert  = 0;
        fgcntt = 0;
        timp   = 0.d0 ;
		
%*----------------------------------------------- Here start minimization               		


 %       call timer(time1)       

      iter = 1;
	  fcnt = 0;
	  gcnt = 0;
      fgcnt = 0; 
      ngama = 0;
%*-----------------------------------------------------------------------

%*                     =====      Beginning the algorithms    =====   
%[fx,grad]= cutest_cons(x);
        [fx,grad]= cutest_obj(x);         %----------
		fcnt=fcnt+1;
		gcnt=gcnt+1;
        fgcnt=fgcnt+2 ;       
        fxmin=fx; 
        fx1 = fx;
        fx2 = fx;
        grad1 = grad;
        grad2 = grad;
		
        ck=fx;
        qk=1.0D+0;
%c        lmda=1  		

%*------------------------------------- Direction
        d = -grad  ;      

	  ginf=norm(grad,Inf);
		
[step,xnew,fxnew,fcnt,fgcnt] = back(ginf,gama,x,d,fx,grad,s,...
                                 fcnt,fgcnt,maxfgcnt);
        [~,gradnew]= cutest_obj(xnew);         %----------
		gcnt=gcnt+1;
        fgcnt=fgcnt+1 ; 
        gradn = norm(grad)^2;
          
           
%* ------------------------------------ Begin Iterations
                             
while true                                                                     
		
		    nck=etak*qk*ck+fxnew;
		    qk=etak*qk+1;
		    ck=nck/qk;
		
        if(fxnew < fxmin) 
            fxmin=fxnew;
        end         
%*-------------------------------------------> gama Barzilai & Borwein

        bs = xnew-x;
        if iter >= 2
            bsp= bs;
            byp= by;
            by = gradnew - grad2;
            bshat= 3./2. *(bs - 1./3. *bsp);
          nuk= 4*fx2 - 1./2. *fx1 - 7./2. *fxnew + 1./2. *(bs+bsp)'*gradnew...
              + bshat'*(4./3. *grad2- 1./3. *grad1);
          zk= by - 1./3. *byp + nuk/(norm(bshat))^2 *bshat;
        else
            by = gradnew - grad;
            zk=by;
        end

          sumgrad = gradnew + grad2;
          sy = bs'*by;
          ss = bs'*bs;
		    yy = by'*by;
		    gs = grad2'*bs;
		    gnews = gradnew'*bs;			
		    sumgs =sumgrad'*bs;

		alpha3 = (2.d0*(fx2-fxnew)+2.d0*gnews)/ss;		
		alpha4 = (6.d0*(fx2-fxnew)+4.d0*gnews+2.d0*gs)/ss;
		alpha11 = (4.d0*(fx2-fxnew)+3.d0*gnews+gs)/ss;
				
	yh = by + ((2.d0*(fx2-fxnew)+sumgs)/ss)*bs;
	yhyh =  yh'*yh;

		alpha5 = yhyh/(6.d0*(fx2-fxnew)+4.d0*gnews+2.d0*gs);		
		alpha6 = yhyh/(2.d0*(fx2-fxnew)+2.d0*gnews);
		alpha12 = yhyh/(4.d0*(fx2-fxnew)+3.d0*gnews+gs);

               phi = 6.d0*(fx2-fxnew) + 3.d0*sumgs;
		       beta = 1.d0 + phi/sy ;
		
		alpha13=beta*sy/ss;
		alpha14=beta*yy/sy;	
        alpha15= (zk'*zk)/(bs'*zk);
			   
       alpha=1.d0/alpha15;			   
		  
      if (alpha < 0.d0)
          alpha=ss/sy;
      end
%      if (alpha .lt. 0.d0)  alpha=sy/yy;
	   	   
% ccc  	   alpha=sy/yy
% ccc        gamab = sy/ss      
% ccc        stepini = 1.d0/gamab				 

       if (alpha < 0.d0) 
	         alpha = lmdamax;
       else 
		     alpha = min(lmdamax,max(lmdamin,alpha));
       end 

		        stepini = alpha	;		
			
%*-------------------------> End gama Barzilai & Borwein
%*------------------------> gama given by Function values
% ccc   if(igama) 
% ccc        if(fxnew-fx+step*gradn < 0.d0)  
% ccc          ngama=ngama+1;
% ccc          vv = (fx - fxnew - step*gradn)/gradn;  
% ccc          eta = delta + vv;
% ccc          step = step + eta;
% ccc        end 
%             
% ccc     gama = 2.d0*(fxnew-fx+step*gradn)/(gradn*step*step)  
% 
% ccc        stepini = 1.d0/gama
% ccc    end   
%*-------------------------------------------> End gama given by Function values      

          d =-gradnew;

%*------------------------------------------------- Backtracking

         step=stepini;
         stepar(iter) = stepini;
         	
% *-----------------------------------------------------------------   
%             
% *---------------------------------------- Updating for a new iteration         

          x = xnew;
          if iter >= 2
            grad1= grad2;
            fx1=fx2;
          end
          grad2 = gradnew;                
          fx2=fxnew;
    
[step,xnew,fxnew,gradnew, fcnt,gcnt,fgcnt ] = backnew( stepini,x,d,...
    gradnew,gama,fcnt,gcnt,fgcnt,ck,s,maxfgcnt )  ;   
%*----------------------------------------------------------------------

 
%* ------------------------------------------------ STOP criterion

        	  ginf=norm(gradnew,Inf);

        if     ginf <= eps 
%             abs(step*ps) <= epsf*abs(fxnew)
           break
        end

        iter=iter+1;
        if(iter > maxiter || fgcnt > maxfgcnt) 
            break
        end

end
%        call timer(time2)
%		        tottime = time2-time1   

                 stepmin = min(stepar);
                 stepmax = max(stepar);
                 
%*---------------------------------------- Write some information

        if 	( iter > maxiter || fgcnt > maxfgcnt) 
            success=0;
	              % WRITE(itfile,140)nexp 
        else
             success=1;
		    %   WRITE(itfile,150)nexp,iter,fcnt,gcnt, tottime 
         end 
      

          itert  = itert  + iter;
          fgcntt = fgcntt + fgcnt ;     
            
%c-------------------------------------- End DO n                             

end