*------------------------------------------------------------Line 1-     
*                      UNCONSTRAINED OPTIMIZATION
*                      ==========================
*                                
*                              Accelerated    
*      Scaled Conjugate Gradient Algorithm BFGS Preconditioned
*                       with Beale-Powell restart 
*      =======================================================
*
*                            Neculai Andrei
*                 Research Institute for Informatics
*              Center for Advanced Modeling and Optimization
*               8-10, Averescu Avenue, Bucharest 1, Romania
*                          E-mail: nandrei@ici.ro
*                             voice: 666.58.74                       
*                           mobile: 0744-777748
*                                                   
*                    Academy of Romanian Scientists
*             54, Splaiul Independentei, Buchrest 5, Romania
*
*                            Main Program
*                           ==============
*
* This is the main program to minimize a train of differentiable 
* functions. The Fortran expression of functions to be minimized
* is given in file "evalfg.for". A particular function is refered 
* by parameter nexp. For example nexp=5 means that the 5-th function
* in the list is considered for minimization.  The algebraic
* expressions of the function and its gradient must be provided
* by the user.
*
* The file "funcname.txt" contains the name of the functions, 
* belonging to this train.    
*
* The initial guess must be introduced in "inipoint" subroutine.
* A particular value for parameter nexp will specify both the
* function to be minimized as well as the initial point. 
*
* The package include a train of 80 functions and their starting guess.
*
* There are two output files where the user can find some
* information about the optimization process. The content of these
* files is self-explained. The user may change the name of these files.   
*
* In case the user want to solve his own problem as an unconstrained
* optimization one he/she have the possibility to place it at the end
* of the EVALFG subroutine. In INIPOINT he/she must introduce the
* initial point where from the computation starts. Besides, he/she
* must introduce the name of the problem into the FUNCNAME.TXT file.
* Pay attention all these must be corelated, i.e. the Fortran 
* expression of the function and its gradients, the initial point and
* the name of the problem must have the same order number (nexp).
*
* The user must change nexpini and nexptot accordingly, and of course
* he/she must change n - the dimension of the problem.        
*
* It is important to have a train of numerical experiments as it is
* organized in this package. You have the possibility to compare
* ASCALCG with some other algorithms, including the Dolan-More'
* profiles plotting.
*
*----------------------------------------------------------------------
*
* Structure of the package:
* =========================
*
* Main Program --------> inipoint 
*                      > gettim
*                      > ASCALCG -------> evalfg
*                                       > LineSearch
*                      > exetim 
*
* Description:
* ============
* inipoint   - Subroutine for initial point specification. 
* gettim     - Subroutine for getting cpu time.   
* exetim     - Subroutine for computing the ellapsed cpu time.
* ASCALCG    - Main Subroutine. 
* evalfg     - Subroutine for function and its gradient evaluation.
* LineSearch - Subroutine for steplength computation. 
*------------------------------------------------------------------
*
*
* Remark
* ------
* This prackege belongs to the FCG Project - "The Fastest Conjugate
* Gradient Algorithm for Nonlinear Unconstrained Optimization".
* It is in close connection to SCALCG package.
*
*
* References
* ----------
* Andrei, N., (2005) SCALCG: A nonlinear scaled conjugate gradient 
*    algorithm for unconstrained optimization. 
*    Technical Report, July 4, 2005 (on the web)
* Andrei, N., (2005) Scaled conjugate gradient algorithms for 
*    unconstrained optimization.
*    Technical Report, March 11, 2005.
* Andrei, N., (2004) A new gradient descent method for unconstrained 
*    optimization.
*    Technical Report, March 12, 2004. 
* Andrei, N., (2006) A scaled nonlinear conjugate gradient algorithm for 
*    unconstrained optimization.
*    OPTIMIZATION, A journal of mathematical programming and operations
*    research. vol. 57, No. 4, August 2008, pp.549-570..
* Andrei, N., (2007) A scaled BFGS preconditioned conjugate gradient algorithm 
*    for unconstrained optimization.
*    Applied Mathematics Letters, 20 (2007), p.645-650.
* Andrei, N., (2006) Scaled memoryless BFGS preconditioned conjugate gradient
*    algorithm for unconstrained optimization.
*    Optimization Methods and Software, vol.22, Number 4, 
*    August 2007, pp.561-571.     
* Andrei, N., (2007) Scaled conjugate gradient algorithms for unconstrained
*    optimization. 
*    Computational Optimization and Applications, vol.38, Number 3, 
*    December 2007, pp.401-416.
* Birgin G.E. and Martinez, J.M., (2001) A spectral conjugate gradient
*    method for unconstrained optimization.
*    Aplied Mathematics and optimization, 43, pp.117-128, 2001.
* Shanno, D.F., (1978) Conjugate gradient methods with inexact 
*    searches.
*    Mathematics of Operations Research, 3, pp.244-256, 1978.
*
*
*
* Date created: March 4, 2008.   
*
* Trace for modifications: 
* March 5, 2008: Implementation of the acceleration scheme
*                in two places immediately after the 
*                line search.
*
*                                              Neculai Andrei, 2008 
*--------------------------------------------------------------------------
*                                           

      real*8 x(500000)
            
      integer iter,fgcnt,lscnt,irstart
      integer itert, irstartt, fgcntt,lscntt 
      integer istopft, istopgt, istopfgt, istopfrt
      integer stoptest                           
      integer n, nexp, nexpini,nexptot, maxiter, igamma
      integer*4 gh,gm,gs,gc, ght,gmt,gst,gct, timpexp    
      integer*2 iyear, imonth, iday

      real*8 epsg,epsf
      real*8 fxnew, gnorm, delta, proc, timp  

      logical tetas, tetaa
      logical accelerat
      character*70, numef, fnumef(200)

* Input file
            
      open(unit=7,file='funcname.txt',status='old')  

* Output files
      
      open(unit=3,file='ascalcg.out',status='unknown')
      open(unit=5,file='ascalcg.rez',status='unknown')


c------------------------ tetas = .true. spectral updating of theta.
*                         tetaa = .true. anticipative updating of theta.

      tetaa  = .false.
      tetas  = .true.                                                     
      
C------------------------ nexptot = total number of problems considered
*                                   in this train of numerical experiments.
      nexpini = 1
      nexptot = 80

c------------------------ stoptest = option parameter for selection of 
*                                    stopping criterion:
*                         if stoptest = 1 then consider the following test:
*                            if(ginf .le. epsg)
*                         if stoptest = 2 then consider the following test:
*                            if(gnorm .le. epsg)
*                         if stoptest = 3 then consider the following test:
*                            if(gnorm .lt. epsg*dmax1(1.d0, dabs(fx)))  
*                         if stoptest = 4 then the stopping criterion is:
*                            if(ginf .le. dmax1(epsg, epsf*ginfz))                        
*                         if stoptest = 5 then the following criteria are used:
*                            if(ginf .le. epsg .OR. dabs(alfa*gtd) .le. epsf*dabs(fx)) 
*                               where:
*                               ginf  = infinite norm of gradient g(xk),
*                               ginfz = infinite norm of gradient g(x0),
*                               gnorm = norm-2 of gradient g(xk).

      stoptest = 1
                       
c------------------------ Activation of the acceleration scheme                       
*                         If accelerat is .true., then the acceleration
*                         scheme is activated. Otherwise the acceleration
*                         scheme is inactive.     
*
      accelerat = .true.
      

c------------------------ Write name of the package and data.

      
      if(accelerat) then 
        write(3,947)                 
947   format(/,2x,'Accelerated SCALCG Algorithm.    Neculai Andrei',/,
     *         2x,'BFGS preconditioned conjugate gradient algorthm',/)

        if(tetas)write(3,48)     
        if(tetaa)write(3,49)  

      else

        write(3,1947)                 
1947  format(/,2x,'SCALCG Algorithm.                Neculai Andrei',/,
     *         2x,'BFGS preconditioned conjugate gradient algorthm',/)

        if(tetas)write(3,48)     
        if(tetaa)write(3,49)	
	end if

        call getdat(iyear, imonth, iday)
        write(3,211) imonth, iday, iyear
211     format(4x,'Date: ---  Month:',i2,'  Day:',i2,'  Year: ',i4,/) 


            
      

*------------------------ Consider the values for tolerances:
*                         epsg: convergence tolerance for gradient.
*                         epsf: convergence tolerance for function values.

      epsg  = 0.000001d0  
      epsf  = 10.d0**(-12)      

*------------------------ Consider the maximum number of iterations and
*                         a value for parameter "delta".
      
      maxiter = 10000        
      delta   = 10.1d0

c------------------------ Read the name of the problems.

              do i=1,80
                read(7,21) numef
                fnumef(i)=numef 
21              format(a70)         
              end do
              

c------------------------ Initializations for total information.
           
        istopft  = 0
        istopgt  = 0  
        istopfgt = 0
        istopfrt = 0
      
*---------------------------------------------  Here Start Experiments   
*                                               ======================
        do nexp = nexpini,nexptot
        
        numef = fnumef(nexp)
        
      if(accelerat) then
        write(3,41) nexp, numef 
        write(*,41) nexp, numef      
41      format(/,2x,i3,2x,'*** ASCALCG Algorithm ***. Function: ',a40)  
      else
        write(3,441) nexp, numef 
        write(*,441) nexp, numef      
441     format(/,2x,i3,2x,'*** SCALCG Algorithm ***. Function: ',a40)  
      
      end if
*
      if(tetas) write(3,43) stoptest
43     format(7x,'Powell criterion for restart.',2x,' Stoptest = ',i1,
     *        2x,'(theta spectral)',/)
      if(tetaa) write(3,44) stoptest        
44     format(7x,'Powell criterion for restart.',2x,' Stoptest = ',i1,
     *        2x,'(theta anticipative)',/) 

*

      if(tetas) write(3,46)
46    format(6x,'n',3x,'iter',4x,'irs',3x,'fgcnt',2x,'lscnt',
     *              3x,'time(c)',7x,'fxnew',16x,'gnorm')
                 
      if(tetaa) write(3,47)
47    format(6x,'n',3x,'iter',4x,'irs',3x,'fgcnt',2x,'lscnt',
     *              3x,'time(c)',7x,'fxnew',16x,'gnorm',
     *              8x,'theta<0')
      
      if(tetas)write(3,48)
48    format(1x,86('-'))      
      if(tetaa)write(3,49)
49    format(1x,93('-')) 

        
c------------------------ Total # of iterations, restart iterations, functions
*                               evaluations, line searches and total cpu time.     
        itert=0
        irstartt=0
        fgcntt=0
        lscntt=0
        timp=0.d0

c--------------------------------------------- Solve the problems
*                                              ==================
*                                           n = number of variables
                        
        do n = 1000, 10000, 1000

          fgcnt=0
          lscnt=0
c------------------------ Select an initial point.
         
         call inipoint(n,x,nexp)

c------------------------ Solve a particular problem.      
*                         ==========================

      call gettim(gh,gm,gs,gc)
       
      call ascalcg(n,x,epsg,epsf,delta,maxiter,stoptest,fxnew,gnorm,
     *             iter,irstart,igamma,fgcnt,lscnt,tetas,tetaa,
     *             nexp,istopf,istopg,istopfg,istopfr,accelerat)
       
      call gettim(ght,gmt,gst,gct)
      call exetim(gh,gm,gs,gc, ght,gmt,gst,gct)
         
c------------------------ Compute elapsed time in centeseconds.

         timpexp = ght*360000 + gmt*6000 + gst*100 + gct 



c------------------------ Write down information in file *.rez.
         
         if(n .eq. 1000) then
           write(5,1010)nexp, n,iter,fgcnt,timpexp, fxnew,gnorm
1010       format(i2,i6,2x,i6,2x,i6,2x,i6,2x,e20.13,2x,e20.13)
         else           
           write(5,101)       n,iter,fgcnt,timpexp, fxnew,gnorm
101        format(2x,i6,2x,i6,2x,i6,2x,i6,2x,e20.13,2x,e20.13)      
         endif
            
            

*------------------------  Write down information in file *.out

      if(tetas) write(3,850)  n, iter, irstart, 
     *                        fgcnt,lscnt,timpexp,
     *                        fxnew, gnorm       
      if(tetas) write(*,1850) n, iter, irstart, 
     *                        fgcnt,lscnt,timpexp,
     *                        fxnew, gnorm       
      if(tetaa) then
        if(igamma .ne. 0) then
          write(3,851)  n, iter, irstart, fgcnt,lscnt,timpexp,
     *                  fxnew,gnorm,igamma       
          write(*,1851) n, iter, irstart, fgcnt,lscnt,timpexp,
     *                  fxnew,gnorm,igamma       
        end if
        if(igamma .eq. 0) then
          write(3,850)  n, iter, irstart, fgcnt,lscnt,timpexp,
     *                  fxnew, gnorm
          write(*,1850) n, iter, irstart, fgcnt,lscnt,timpexp,
     *                  fxnew, gnorm  
        end if
      end if   
      
850   format(1x,i6,1x,i6,1x,i6,2x,i6,1x,i6,1x,i7,1x,
     *                         e20.13,1x,e20.13)  
1850  format(1x,i6,1x,i6,1x,i6,2x,i6,1x,i6,1x,i7,1x,
     *                         e15.8,1x,e15.8)      
851   format(1x,i6,1x,i6,1x,i6,2x,i6,1x,i6,1x,i7,1x,
     *                         e20.13,1x,e20.13,1x,i5)
1851  format(1x,i6,1x,i6,1x,i6,2x,i6,1x,i6,1x,i7,1x,
     *                         e15.8,1x,e15.8,1x,i5)

*------------------------ End Write
                                                              
          itert    = itert    + iter
          irstartt = irstartt + irstart
          fgcntt   = fgcntt   + fgcnt
          lscntt   = lscntt   + lscnt
          timp     = timp     + float(timpexp) 

          istopft  = istopft  + istopf
          istopgt  = istopgt  + istopg            
          istopfgt = istopfgt + istopfg
          istopfrt = istopfrt + istopfr             

c-------------------------------------------------- End DO n         
          end do 

          if(itert .ne. 0) then
            proc=float(irstartt)*100.d0/float(itert)                                                     
          else
            proc = 0.d0  
          end if  
                                                     
       if(tetas)write(3,48)        
       if(tetaa)write(3,49)
       write(3,55) itert, irstartt, fgcntt, lscntt,timp/100.d0, proc
55     format(1x,'TOTAL',1x,i7,1x,i6,1x,i7,1x,i6,1x,f7.2,'(seconds)',
     *                   5x,'proc= ',f5.2,'%',/)   
     
       write(*,755) timp/100.d0
755    format(1x,'Total cpu time = ',f7.2)
       
c---------------------------------------------------------- End do nexp                         
      end do    


      write(5,18)
18    format(2x,'     n',2x,'  iter',2x,' fgcnt',2x,'time(c)')  

      if(stoptest .eq. 5) then      
        write(3,1751) istopft, istopgt, istopfgt, istopfrt
        write(5,1751) istopft, istopgt, istopfgt, istopfrt                   
1751    format(/,1x,'SCALCG: istop (f)  total =',i7,
     *                   4x,'istop (g)  total =',i7,/,
     *                   9x,'istop (fg) total =',i7,
     *                   4x,'istop (fr) total =',i7) 
      end if

      write(5,947)  
      if(tetas) write(5,143) stoptest
143     format(2x,'Powell criterion for restart.',2x,' Stoptest = ',i1,
     *        2x,'(theta spectral)',/)
      if(tetaa) write(5,144) stoptest        
144     format(2x,'Powell criterion for restart.',2x,' Stoptest = ',i1,
     *        2x,'(theta anticipative)',/) 
     
      stop
      end
*----------------------------------------------- END of Main Program
                                                  





*                   *****  ASCALCG  *****  ASCALCG  *****     
*                   =====================================
*
*-------------------------------------------------------------------
*                           Subroutine ASCALCG                                    
*                           ==================
*                                        
*                              Accelerated
*                 Scaled Conjugate Gradient Algorithms for
*                       Unconstrained Optimization
*
*
*                            Neculai Andrei
*                 Research Institute for Informatics
*              Center for Advanced Modeling and Optimization
*               8-10, Averescu Avenue, Bucharest 1, Romania
*                          E-mail: nandrei@ici.ro
*                             voice: 666.58.70                       
*                                 and           
*                    Academy of Romanian Scientists
*             54, Splaiul Independentei, Buchrest 5, Romania                       
*
* /-----------------------------------------------------------------\
* | ASCALCG is a subroutine dedicated to compute the minimizer of a |
* | differentiable function with a large number of variables.       |
* | ASCALCG implements a nonlinear scaled conjugate gradient        |
* | algorithm.                                                      |
* | This subroutine is accompanied by subroutine "LineSearch" which |
* | implements the Wolfe line search. Both these subroutines belong |
* | to ASCALCG package.                                             |
* |                                                                 |
* | The user must provide a subroutine to evaluate the function     |
* | value and its gradient in a point. The name of this subroutine  |
* | is EVALFG. On the web there is a file "func.for" where over 70  |
* | Fortran expression of some test functions are included.         |
* | Additionally, the algebraic expression of these functions is    |
* | given in file "testuo.pdf".                                     |
* |                                                                 |
* | There are some facilities for the user to specify:              |
* | 1) The termination criterion.                                   |
* | 2) Convergence tolerance for gradient.                          |
* | 3) Convergence tolerance for function value.                    |
* | 4) Maximum number of iterations in SCALCG subroutine.           |
* | 5) Maximum number of iterations in LineSearch subroutine.       |
* | 6) Level of printing.                                           |
* |                                                                 |
* |-----------------------------------------------------------------|      
* |                                                                 |
* | The calling sequence of ASCALCG is:                             |
* |                                                                 |
* |   call ascalcg(n,x,epsg,epsf,delta,maxiter,stoptest,fx,gnorm,   |
* |                iter,irstart,igamma,fgcnt,lscnt,tetas,tetaa,     |
* |                nexp, istopf,istopg,istopfg,istopfr,accelerat)   |
* |                                                                 |
* |                                                                 |
* |Input parameters:                                                |
* |=================                                                |
* |n          (integer) number of variables.                        |
* |x          (double)  starting guess, length n. On output         |
* |                     contains the solution.                      |
* |epsg       (double)  convergence tolerance for gradient.         |
* |epsf       (double)  convergence tolerance for function.         |
* |delta      (double)  parameter used in anticipative selection    |
* |                     the scaling factor of gradient. It has a    |
* |                     value comparable with the function's        |
* |                     value in the current point.                 |
* |maxiter    (integer) maximum number of iterations.               |
* |stoptest = option parameter for selection of                     | 
* |            stopping criterion:                                  |
* |            if stoptest = 1 then consider the following test:    | 
* |               if(ginf .le. epsg)                                |
* |            if stoptest = 2 then consider the following test:    |
* |               if(gnorm .le. epsg)                               |
* |            if stoptest = 3 then consider the following test:    |
* |               if(gnorm .lt. epsg*dmax1(1.d0, dabs(fx)))         |
* |            if stoptest = 4 then the stopping criterion is:      |
* |               if(ginf .le. dmax1(epsg, epsf*ginfz))             |
* |            if stoptest = 5 then the following criteria are used:|
* |               if(ginf .le. epsg .OR.                            |
* |                  dabs(alfa*gtd) .le. epsf*dabs(fx))             |
* |               where:                                            |
* |               ginf  = infinite norm of gradient g(xk),          | 
* |               ginfz = infinite norm of gradient g(x0),          |
* |               gnorm = norm-2 of gradient g(xk).                 |
* |tetas      (logical) if tetas = .true. then the spectral formula |
* |                     is used by the algorithm.                   |
* |tetaa      (logical) if tetaa = .true. then the anticipative     |
* |                     formula is used by the algorithm.           |
* |nexp       (integer) parameter specifying the number of the      |
* |                     problem considered in a train of experiments|
* |accelerat  (logical) parameter for activation of the acceleration|
* |                     scheme. If accelerat = .true., then the     |
* |                     acceleration scheme is activated.           |
* |                                                                 |
* |                                                                 |
* |Output parameters:                                               |
* |==================                                               |
* |fx         (double)  function value in final (optimal) point.    |
* |gnorm      (double)  norm-2 of gradient in final point.          |
* |iter       (integer) number of iterations to get the final point.|
* |irstart    (integer) number of restart iterations.               |
* |igamma     (integer) number of gamma negative in anticipative    |
* |                     formula.                                    |
* |fgcnt      (integer) number of function evaluations.             |
* |lscnt      (integer) number of line searches.                    |
* |istopf     (integer) = 1 if "dabs(alfa*gtd) .le. epsf*dabs(fx)", |
* |                     = 0 othewise.                               |
* |istopg     (integer) = 1 if "infinit||grad|| .le. epsg",         |
* |                     = 0 otherwise.                              |
* |istopfg    (integer) = 1 if "dabs(alfa*gtd) .le. epsf*dabs(fx)"  |
* |                            .AND.                                |
* |                            "infinit||grad|| .le. epsg",         |
* |                     = 0 otherwise.                              |
* |istopfr    (integer) = 1 if "dabs(fxnew-fx)/(1.d0+dabs(fx))      |
* |                             .le. epsf" in the segment of code   |
* |                         for theta computed in an anticipative   |
* |                         manner.                                 |
* |-----------------------------------------------------------------|
* |                                                                 |
* |                                                                 |
* |Calling subroutines:                                             |
* |====================                                             | 
* |Subroutine ASCALCG is calling two subroutines:                   |
* |EVALFG     an user subroutine, and                               |
* |LINESEARCH a package subroutine.                                 |
* |                                                                 |   
* |The user must supply a subroutine with the function and its      |
* |gradient:                                                        |
* |  call evalfg(n,x,fx,grad,nexp)                                  |
* |where:                                                           |
* |  n    (integer)  number of variables.                           |
* |  x    (double)   the current iteration.                         |
* |  fx   (double)   function value in point x.                     |
* |  grad (double)   array with gradient of function in point x.    |
* |  nexp (integer)  parameter specifying the number of the         |
* |                  problem considered in a train of experiments.  |
* |                                                                 |
* |                                                                 |
* |The Wolfe line search is implemented in the subroutine           |
* |LINESEARCH, which belongs to the package.                        |
* |The calling sequence of this subroutine is as follows:           |            
* | call LineSearch (n,x,f,d,gtd,dnorm,alpha,xnew,fnew,gnew,        |
* |                  fgcnt,lscnt,lsflag, nexp)                      |
* |where:                                                           |
* |  n      (integer)  number of variables.                         |
* |  x      (double)   the current iteration.                       |
* |  f      (double)   function value in current point.             |
* |  d      (double)   array with search direction.                 |
* |  gtd    (double)   scalar: grad'*d.                             |   
* |  dnorm  (double)   2 norm of d.                                 |
* |  alpha  (double)   step length (given by the LineSearch).       |
* |  xnew   (double)   array with the new estimation of variables.  |
* |  fnew   (double)   function value in xnew.                      |
* |  gnew   (double)   array with gradient in xnew.                 |
* |  fgcnt  (integer)  number of function evaluations.              |
* |  lscnt  (integer)  number of line searches.                     |
* |  lsflag (integer)  parameter for abnormal Line Search           |
* |                    Termination. If the # of iterations in       |
* |                    LineSearch is greater than 20 then lsflag=1  |
* |  nexp   (integer)  parameter specifying the number of the       |
* |                    problem considered in a train of experiments |
* |                                                                 |
* |Subroutine LINESEARCH is the same as that used in CONMIN package |
* |by Shanno and Phua, and in SCG package by Birgin and Martinez.   |
* |                                                                 |
* |-----------------------------------------------------------------|
* |                                     Neculai Andrei, 2005/2008   |
* \-----------------------------------------------------------------/
           
      subroutine ascalcg(n,x,epsg,epsf,delta,maxiter,stoptest,fx,
     *             gnorm,iter,irstart,igamma,fgcnt,lscnt,tetas,tetaa,
     *             nexp, istopf,istopg,istopfg,istopfr, accelerat)  
*
      parameter(ia=500000) 
      
      integer n, stoptest, fgcnt,lscnt,lsflag  
      integer istopf,istopg,istopfg,istopfr        
      integer iter, irstart, igamma    
      integer maxiter, nexp
      real*8  epsf, epsg
      real*8  alfa, alfan              
      real*8  gnorm, dnorm, dnormpre, ginf, ginfz
      real*8  fx, fxnew, gtg, g1tg, gtd
      real*8  theta, thetas, gtdr, dtd 
      real*8  gts,gty,yts,yty,yks,yky,gtw,ytw, coef, sts,sty
      real*8  gts1,yts1 , an,bn
      real*8  gamma, eta, delta 
      real*8 x(ia), xnew(ia), grad(ia), gradnew(ia), d(ia)
      real*8 s(ia),y(ia), ss(ia), ys(ia), v(ia),w(ia)
      logical tetas, tetaa     
      logical accelerat
      
* Step 1. Initialization.  Initial step
                
      iter=0   
      irstart=0           
      igamma=0

          istopf  = 0
          istopg  = 0           
          istopfg = 0
          istopfr = 0
                
      call evalfg(n,x,fx,grad,nexp) 
      fgcnt=fgcnt+1

      gtg   = 0.d0
      ginfz = 0.d0
      do i=1,n
        d(i) = -grad(i)
        gtg = gtg + grad(i)**2
        ginfz = dmax1(ginfz, dabs(grad(i)))
      end do                                            
      
*         ginfz = infinite norm of grad(x0). It is used in stopping
*                 test of iterations, for stoptest .eq. 4.

      gtd = -gtg
      gnorm = sqrt(gtg)
      dnorm = gnorm
      gtdr = gtd
      dtd  = gtg

      
            
* Step 2. Line search from initial point. Determine the step length alfa.
*         Initialize alfa and then with this value compute the step length.
*         Compute xnew = x + alfa*d,  fnew = f(xnew),  gradnew = grad(xnew).

      if(gnorm .ne. 0.d0) alfa=1.d0/gnorm          

      call linesearch(n,x,fx,d,gtd,dnorm,alfa, xnew,fxnew,gradnew,
     *                fgcnt,lscnt,lsflag, nexp)       

c      if(lsflag .ne. 0) then
c        write(*,*) 'Warning ! Abnormal Line Search Termination.'      
c      end if



* Step 3. Compute s and y. Prepare for test of continuation. 
*                          ginf = infinite norm of gradient
*                          gtg = gradnew'*gradnew (scalar product)
*                          gnorm = norm 2 of gradient
        
      gtg  = 0.d0
      ginf = 0.d0
      do i=1,n
        s(i) = xnew(i) - x(i)
        y(i) = gradnew(i) - grad(i) 
        gtg  = gtg + gradnew(i)*gradnew(i)
        ginf = dmax1(ginf, dabs(gradnew(i)))
      end do  
      
      gnorm = sqrt(gtg)               



               
* Step 4. Stopping test for initial point. This test is used only here. 
*         Inhibited.

c      if(ginf .lt. epsg .or. 
c     *   dabs(fxnew-fx)/(1.d0+dabs(fx)) .le. epsf) go to 999 


                                                                        
                                                                        
* Step 5. Continue iterations and test for maximum number of iterations.

              iter = iter + 1
              if(iter .gt. maxiter) then
                go to 1000
              end if          

            

111   continue
*--------------------------- Here a restart iteration,
*                            i.e. the Beale-Powell restart criterion has been satisfied.




* Step 6. Compute theta in a spectral or an anticipative manner. 
*
*                                           -----------------------
*------------------------------------------ | 1) Spectral Formula |
*                                           -----------------------
      
      if(tetas) then   
      
        sts = 0.d0
        sty = 0.d0 
        do i=1,n
          sts = sts + s(i)*s(i)
          sty = sty + s(i)*y(i)  
        end do                
c                        
        if(sty .le. 0.d0) then
          theta = 1.0d+10
        else
c          theta = dmin1(1.0d+10, dmax1(1.0d-10, sts/sty))
          theta = sts/sty
        end if 
 
      end if           

                          
*                                            ---------------------------
*------------------------------------------- | 2) Anticipative Formula | 
*                                            ---------------------------
      
      if(tetaa) then
                                                            
c          if( dabs(fxnew-fx)/(1.d0+dabs(fx)) .le. epsf) then
c            istopfr=istopfr + 1
c            go to 999                                                                 
c          end if  
                                                            
          gamma = 2.d0*(fxnew - fx - alfa*gtdr)/(dtd*alfa*alfa)   
        
        if(gamma .le. 0.d0) then       
          igamma = igamma+1
          eta   = (fx-fxnew+alfa*gtdr+delta)/gtdr  
          alfan = alfa - eta
          gamma = 2.d0*(fxnew - fx - alfan*gtdr)/(dtd*alfan*alfan)     
        end if 
        
        theta = 1.d0/gamma
      end if

*------------------------------------------- End theta computation

* Save these values.
                          
      fx=fxnew
      dnormpre = dnorm




                          
* Step 7. Compute direction d.  
*                Here 4 scalar products for restart d(k+1). 
*                Firstly update x and grad.
*                Line search preparations for this new direction.
*                Compute the scalar products: gtd and dtd.

      gts=0.d0
      gty=0.d0
      yts=0.d0
      yty=0.d0                             
      do i=1,n
          x(i) = xnew(i)
          grad(i) = gradnew(i)
        gts = gts + grad(i)*s(i)
        gty = gty + grad(i)*y(i)
        yts = yts + y(i)*s(i)
        yty = yty + y(i)*y(i)
      end do
                  
      if(yts .ne. 0.d0) then
        coef = 1.d0 + theta*yty/yts                  
        gtd = 0.d0
        dtd = 0.d0
        do i=1,n
          d(i) = -theta * grad(i) 
     *           +theta * (gts/yts) * y(i)  
     *           -( coef*gts/yts - theta*gty/yts ) * s(i)      
     
          gtd = gtd + grad(i)*d(i)
          dtd = dtd + d(i)*d(i)
        end do
      end if
      
      if(yts .eq. 0.d0) then   
        gtd = 0.d0
        dtd = 0.d0
        do i=1,n
          d(i) = -theta * grad(i)
          gtd = gtd + grad(i)*d(i)
          dtd = dtd + d(i)*d(i)
        end do
      end if      
      
      
      dnorm = sqrt(dtd)
      gtdr = gtd




* Step 8. Line search.
*         Initialize alfa and then with this value compute the step length.
*         Compute xnew = x + alfa*d,  fnew = f(xnew),  gradnew = grad(xnew).

      alfa = alfa * (dnormpre/dnorm)

      call linesearch(n,x,fx,d,gtd,dnorm,alfa, xnew,fxnew,gradnew,
     *                fgcnt,lscnt,lsflag, nexp)       

c      if(lsflag .ne. 0) then
c        write(*,*) 'Warning ! Abnormal Line Search Termination.'      
c      end if      


                    
                    
c Step 9. Acceleration scheme (I)

      if(accelerat) then      
        an=0.d0
        bn=0.d0
        do i=1,n
          an = an + grad(i)*d(i)
          bn = bn + (grad(i)-gradnew(i))*d(i)
        end do
      
        if(bn .ne. 0.d0) then
          do i=1,n
            xnew(i) = x(i) + (an/bn)*alfa*d(i)
          end do
          call evalfg(n,xnew,fxnew,gradnew, nexp) 
          fgcnt=fgcnt+1
        end if      
      end if


      
      
* Step 10. Store theta
*
                   
      thetas=theta
        


        
* Step 11. Compute s and y. Store these vectors for a normal step.
*          Prepare for stopping test and for Powell restart test.
*          Compute ginf for its usage in stopping test.
*          Compute gtg for its usage in Powell restart test.

      gtg  = 0.d0
      ginf = 0.d0                         
      do i=1,n
        s(i) = xnew(i) - x(i)
        y(i) = gradnew(i) - grad(i)      
        ss(i) = s(i)
        ys(i) = y(i)  
        gtg = gtg + gradnew(i)*gradnew(i)
        ginf = dmax1(ginf, dabs(gradnew(i)))
      end do
      
      gnorm = sqrt(gtg)




*                           
* Step 12. Test of Stopping the iterations  ---------------------- (I) 
*
      if(stoptest .eq. 1) then
        if(ginf .le. epsg) go to 999
      end if  

      if(stoptest .eq. 2) then
        if(gnorm .le. epsg) go to 999
      end if 

      if(stoptest .eq. 3) then
        if(gnorm .le. epsg * dmax1(1.d0, dabs(fxnew))) go to 999
      end if  

      if(stoptest .eq. 4) then
        if(ginf .le. dmax1(epsg, epsf*ginfz)) go to 999
      end if
      
      if(stoptest .eq. 5) then
        if(ginf .le. epsg .or.
     *    dabs(alfa*gtd) .le. epsf*dabs(fxnew)) go to 999 
      end if



      
* Step 13. Continue iterations.

222   continue
 
              iter = iter + 1
              if(iter .gt. maxiter) then
                go to 1000
              end if                      



*      
* Step 14. Test for restart.      **** Powell criterion. ****             
*                                 ===========================
        g1tg=0.d0
        do i=1,n
          g1tg = g1tg + gradnew(i)*grad(i)            
        end do
      
        if(dabs(g1tg) .ge. 0.2d0 * gtg) then 
          irstart = irstart + 1
          go to 111
        end if  
                                                                            


                                                                            
* Step 15.       
*--------------------------------- Here a Standard iteration.
*                                  The Beale-Powell restart criterion is not satisfied.
*                                  Compute d(k+1) using the BFGS philosophy.                   

      fx=fxnew           
      dnormpre=dnorm           


      

* Step 16. Compute v, w, d(k+1)
*                                  Here 6 scalar products for v and w 
*                                       2 scalar products for d(k+1): gts1 and yts1                   
      gts=0.d0
      gty=0.d0
      yts=0.d0
      yty=0.d0
      yks=0.d0
      yky=0.d0       
        gts1=0.d0
        yts1=0.d0
      do i=1,n 
          x(i)    = xnew(i)
          grad(i) = gradnew(i)      
        gts = gts + grad(i)*ss(i)
        gty = gty + grad(i)*ys(i)
        yts = yts + ys(i)*ss(i)
        yty = yty + ys(i)*ys(i)
        yks = yks + y(i) *ss(i)
        yky = yky + y(i) *ys(i) 
          gts1 = gts1 + grad(i)*s(i)
          yts1 = yts1 + y(i)*s(i)
      end do       
      
      coef = 1.d0 + thetas*yty/yts
      
      gtw=0.d0 
      ytw=0.d0
      
      do i=1,n
        v(i) = thetas*grad(i) 
     *        -thetas*(gts/yts)*ys(i) 
     *        +(coef*gts/yts - thetas*gty/yts)*ss(i)
     
        w(i) = thetas*y(i)
     *        -thetas*(yks/yts)*ys(i) 
     *        +(coef*yks/yts - thetas*yky/yts)*ss(i)

        gtw = gtw + grad(i)*w(i) 
        ytw = ytw + y(i)*w(i)

      end do
*
*                                  Here compute standard d(k+1)
*                                  Prepare for line search for this direction:
*                                  compute gtd and dtd.
      gtd = 0.d0
      dtd = 0.d0        
      do i=1,n
        d(i) = -v(i)
     *         +(gts1*w(i) + gtw*s(i))/yts1                 
     *         -(1.d0+ytw/yts1)*(gts1/yts1)*s(i)
     
        gtd = gtd + grad(i)*d(i)
        dtd = dtd + d(i)*d(i)
      end do

      dnorm = sqrt(dtd)
      gtdr =gtd





* Step 17. Line search.
*          Initialize alfa and then with this value compute the step length.
*          Compute xnew = x + alfa*d,  fnew = f(xnew),  gradnew = grad(xnew).


      alfa = alfa * (dnormpre/dnorm)

      call linesearch(n,x,fx,d,gtd,dnorm,alfa, xnew,fxnew,gradnew,
     *                fgcnt,lscnt,lsflag, nexp)       

c      if(lsflag .ne. 0) then
c        write(*,*) 'Warning ! Abnormal Line Search Termination.'      
c      end if      





c Step 18. Acceleration scheme  (II)

      if(accelerat) then
        an=0.d0
        bn=0.d0
        do i=1,n
          an = an + grad(i)*d(i)
          bn = bn + (grad(i)-gradnew(i))*d(i)
        end do
      
        if(bn .ne. 0.d0) then
          do i=1,n
            xnew(i) = x(i) + (an/bn)*alfa*d(i)
          end do
          call evalfg(n,xnew,fxnew,gradnew, nexp)
          fgcnt=fgcnt+1
        end if      
      end if




* Step 19. Compute s and y.
*          Prepare for stopping test and for Powell restart test.
*          Compute ginf for its usage in stopping test.
*          Compute gtg for its usage in Powell restart test.

      gtg  = 0.d0
      ginf = 0.d0
      do i=1,n
        s(i) = xnew(i) - x(i)
        y(i) = gradnew(i) - grad(i)      
        gtg  = gtg + gradnew(i)*gradnew(i) 
        ginf = dmax1(ginf, dabs(gradnew(i)))
      end do                                
      
      gnorm = sqrt(gtg)
                           



*                                 
* Step 20. Test of Stopping the iterations  ---------------------- (II)
*

      if(stoptest .eq. 1) then
        if(ginf .le. epsg) go to 999
      end if  

      if(stoptest .eq. 2) then
        if(gnorm .le. epsg) go to 999
      end if 

      if(stoptest .eq. 3) then
        if(gnorm .le. epsg * dmax1(1.d0, dabs(fxnew))) go to 999
      end if  

      if(stoptest .eq. 4) then
        if(ginf .le. dmax1(epsg, epsf*ginfz)) go to 999
      end if
      
      if(stoptest .eq. 5) then
        if(ginf .le. epsg .or.
     *    dabs(alfa*gtd) .le. epsf*dabs(fxnew)) go to 999 
      end if
      



* Step 21. Continue iterations.
*

      go to 222      



      
* Step 22. Here criteria of STOPping the iterations are satisfied.

999   continue       

         
* Test which criterion stopped the iterations.

      if(stoptest .eq. 5) then                        
        if(ginf .le. epsg) istopg = 1                        
        if(dabs(alfa*gtd) .le. epsf*dabs(fx)) istopf = 1
        if(ginf .lt. epsg .and.
     *    (dabs(alfa*gtd) .le. epsf*dabs(fx))) istopfg = 1   
      end if       
        
c------------------------ Write down the solution, if the case.
                        
c      write(9,802)(i,x(i),i=1,n)
c802   format(10x,i5,4x,e20.13)       


1000  continue

      return
      end
*-------------------------------------------------- End ASCALCG            
      
      
      
      


      
* -------------------------------------------------------------------
* |                     Subroutine LineSearch                       |
* |                     =====================                       |
* |                                                                 |
* |This subroutine implements the Wolfe line search conditions.     |
* |LineSearch subroutine belongs to the package.                    |
* |                                                                 |
* |The calling sequence of this subroutine is as follows:           |            
* |                                                                 |
* | call LineSearch (n,x,f,d,gtd,dnorm,alpha,xnew,fnew,gnew,        |
* |                  fgcnt,lscnt,lsflag, nexp)                      |
* |                                                                 |
* |Input Parameters:                                                |
* |=================                                                |
* |  n      (integer)  number of variables.                         |
* |  x      (double)   the current iteration.                       |
* |  f      (double)   function value in current point.             |
* |  d      (double)   array with search direction.                 |
* |  gtd    (double)   scalar: grad'*d.                             |   
* |  dnorm  (double)   2 norm of d.                                 |
* |  nexp   (integer)  parameter specifying the number of the       |
* |                    problem considered in a train of experiments |
* |Output Parameters:                                               |
* |==================                                               |
* |  alpha  (double)   step length (given by the LineSearch).       |
* |  xnew   (double)   array with the new estimation of variables.  |
* |  fnew   (double)   function value in xnew.                      |
* |  gnew   (double)   array with gradient in xnew.                 |
* |  fgcnt  (integer)  number of function evaluations.              |
* |  lscnt  (integer)  number of line searches.                     |
* |  lsflag (integer)  parameter for abnormal Line Search           |
* |                    Termination. If the # of iterations in       |
* |                    LineSearch is greater than 20 then lsflag=1  |
* |                                                                 |
* |-----------------------------------------------------------------|


      subroutine LineSearch (n,x,f,d,gtd,dnorm,alpha,xnew,fnew,gnew,
     +                       fgcnt,lscnt,lsflag, nexp)

C This is the one-dimensional line search used in CONMIN by Shanno and Phua.
C Please see the subroutine LineSearch in SCG package by Birgin and Martinez.

      integer n,fgcnt,lscnt,lsflag
      double precision f,gtd,dnorm,alpha,fnew
      double precision x(n),d(n),xnew(n),gnew(n)
      integer i,lsiter
      double precision alphap,alphatemp,fp,dp,gtdnew,a,b

      lsflag = 0

* Maximum number of LineSearch is max$ls (now=20)
* Remark: Since ASCALCG implements the acceleration technique (if accelerated = .true.)
* it follows that we may reduce the maximum number of linesearch, let say 2 or 3.
* Therefore, when accelerated = .true., then put max$ls = 2 (or 3). 

      max$ls=20
      
      alphap = 0.0d0
      fp     = f
      dp     = gtd

      do i = 1,n
          xnew(i) = x(i) + alpha * d(i)
      end do
c*1
      call evalfg(n,xnew,fnew,gnew, nexp)
      fgcnt = fgcnt + 1

      gtdnew = 0.0d0
      do i = 1,n
          gtdnew = gtdnew + gnew(i) * d(i)
      end do

      lsiter = 0

* Test whether the Wolfe line seach conditions have been met.

 10   if ( alpha * dnorm .gt. 1.0d-30 .and. lsiter .lt. max$ls .and.
     +     .not. ( gtdnew .eq. 0.0d0 .and. fnew .lt. f ) .and.
     +     ( ( fnew .gt. f + 1.0d-04 * alpha * gtd .or. 
     +     abs( gtdnew / gtd ) .gt. 0.9d0 ) .or. ( lsiter .eq. 0 .and.
     +     abs( gtdnew / gtd ) .gt. 0.5d0 ) ) ) then

* Thest whether the new point has a negative slope and a higher function
* value than that corresponding to alpha=0. In this case, the search has 
* passed through a local max and is heading for a distant local minimum.
* Reduce alpha, compute the new corresponding point, its function value 
* and gradient, as well as gtdnew.
* Repeat this test until a good point has been found.

 20       if ( alpha * dnorm .gt. 1.0d-30 .and. fnew .gt. f .and.
     +         gtdnew .lt. 0.0d0 ) then

              alpha = alpha / 3.0d0

              do i = 1,n
                  xnew(i) = x(i) + alpha * d(i)
              end do
c*2
              call evalfg(n,xnew,fnew,gnew, nexp)
              fgcnt = fgcnt + 1

              gtdnew = 0.0d0
              do i = 1,n
                  gtdnew = gtdnew + gnew(i) * d(i)
              end do

              alphap = 0.0d0
              fp     = f
              dp     = gtd

              go to 20

          end if

* Cubic interpolation to find a new trial point corresponding to alphatemp.

          a = dp + gtdnew - 3.0d0 * ( fp - fnew ) / ( alphap - alpha )
          b = a ** 2 - dp * gtdnew

          if ( b .gt. 0.0d0 ) then
              b = sqrt( b )
          else
              b = 0.0d0
          end if

          alphatemp = alpha - ( alpha - alphap ) * ( gtdnew + b - a ) /
     +                ( gtdnew - dp + 2.0d0 * b )

* Test whether the line minimum has been bracketed.

          if ( gtdnew / dp .le. 0.0d0 ) then

* Here the minimum has been bracketed.
* Test whether the trial point lies sufficiently within the bracketed 
* interval. If it does not, choose alphatemp as the midpoint of the
* interval.

              if ( 0.99d0 * max( alpha, alphap ) .lt. alphatemp .or.
     +            alphatemp .lt. 1.01d0 * min( alpha, alphap ) ) then
                  alphatemp = ( alpha + alphap ) / 2.0d0
              end if

          else

* Here the minimum has not been bracketed.
* The trial point is too small, double the largest prior point.

              if ( gtdnew .lt. 0.0d0 .and. 
     +            alphatemp .lt. 1.01d0 * max( alpha, alphap ) ) then
                  alphatemp = 2.0d0 * max( alpha, alphap )
              end if

* The trial point is too large, halve the smallest prior point.

              if ( ( gtdnew .gt. 0.0d0 .and.
     +            alphatemp .gt. 0.99d0 * min( alpha, alphap ) ) .or.
     +            alphatemp .lt. 0.0d0 ) then
                  alphatemp = min( alpha, alphap ) / 2.0d0
              end if

          end if

* Save and continue the search.

          alphap = alpha
          fp     = fnew
          dp     = gtdnew

          alpha = alphatemp

          do i = 1,n
              xnew(i) = x(i) + alpha * d(i)
          end do
c*3
          call evalfg(n,xnew,fnew,gnew, nexp)
          fgcnt = fgcnt + 1

          gtdnew = 0.0d0
          do i = 1,n
              gtdnew = gtdnew + gnew(i) * d(i)
          end do

          lsiter = lsiter + 1

      go to 10                                                       
      
c-
      end if

      if ( lsiter .ge. max$ls ) then
          lsflag = 1
      end if

      if ( lsiter .ne. 0 ) then
          lscnt = lscnt + 1
      end if

      return

      end
*----------------------------------------------- End of LineSearch
                                
                                
c         End of ASCALCG package: (ASCALCG & LineSearch) 





                                
                                
*----------------------------------------------------------------
*  Date created       : May 30, 1995
*  Date last modified : May 30, 1995
*
*  Subroutine for execution time computation.
*
*----------------------------------------------------------------
*
	subroutine exetim(tih,tim,tis,tic, tfh,tfm,tfs,tfc)
*
	  integer*4 tih,tim,tis,tic
	  integer*4 tfh,tfm,tfs,tfc
*
	  integer*4 ti,tf
	  integer*4 ch,cm,cs
	  data ch,cm,cs/360000,6000,100/
*
	  ti=tih*ch+tim*cm+tis*cs+tic
	  tf=tfh*ch+tfm*cm+tfs*cs+tfc
	  tf=tf-ti
	  tfh=tf/ch
	  tf=tf-tfh*ch
	  tfm=tf/cm
	  tf=tf-tfm*cm
	  tfs=tf/cs
	  tfc=tf-tfs*cs
*
	  return
	end
*-------------------------------------------------- End of EXETIM











*
*
* .-----------------------------------------------------------------.
* |                     Subroutine inipoint              80 problems|
* |                     ===================                         |
* |                                                                 |
* |Subroutine for initial point specification.                      |
* |This is a user subroutine:                                       |
* |                                                                 |
* | The calling sequence is:                                        |
* |                                                                 |
* |  call inipoint(n,x,nexp)                                        |
* |where:                                                           |
* |  n    (integer) the number of variables,                        |
* |  x    (double)  array with the initial point.                   |
* |  nexp (integer) parameter specifying the number of the          |
* |                 problem considered in a train of experiments.   |  
* |                                                                 |
* |                                                  Neculai Andrei |
* .-----------------------------------------------------------------.


******************************************************************
      subroutine inipoint(n,x, nexp)

C     This subroutine computes the initial point

      real*8 x(n)
      
      go to ( 1, 2, 3, 4, 5, 6, 7, 8, 9,10,
     *       11,12,13,14,15,16,17,18,19,20,
     *       21,22,23,24,25,26,27,28,29,30,
     *       31,32,33,34,35,36,37,38,39,40,
     *       41,42,43,44,45,46,47,48,49,50,
     *       51,52,53,54,55,56,57,58,59,60,
     *       61,62,63,64,65,66,67,68,69,70,
     *       71,72,73,74,75,76,77,78,79,80) nexp

1     continue
c                             Freudenstein & Roth - FREUROTH (CUTE)
        i=1
91      x(i)  =   0.5d0
        x(i+1)=  -2.d0
        i=i+2
        if(i.le.n) go to 91      
      return
      
2     continue
c                             Trigonometric        
      do i=1,n
        x(i) = 0.2d0 
      end do   
      return
      
3     continue
c                             Extended Rosenbrock  SROSENBR (CUTE)
        i=1
93      x(i)  =  -1.2d0
        x(i+1)=   1.d0
        i=i+2
        if(i.le.n) go to 93      
      return



4     continue
c                             Extended White & Holst
        i=1
94      x(i)  =  -1.2d0
        x(i+1)=   1.d0
        i=i+2
        if(i.le.n) go to 94      
      return                


5     continue
c                             Extended Beale
        i=1
95      x(i)  =   1.d0
        x(i+1)=   0.8d0
        i=i+2
        if(i.le.n) go to 95      
      return      
      
6     continue
c                             Penalty
      do i=1,n  
       x(i) = float(i)        
      end do   
      return
      
7     continue
c                             Perturbed Quadratic
      do i=1,n
        x(i) = 1.d0 
      end do   
      return      
      
8     continue
c                             Raydan 1
      do i=1,n
        x(i) = 1.d0 
      end do   
      return     
      
9     continue
c                             Raydan 2
      do i=1,n
        x(i) = 1.d0 
      end do   
      return       


10    continue
c                             TR-SUMM
      do i=1,n
        x(i) = 1.d0 
      end do   
      return                    

      
11    continue
c                             Diagonal 1
      do i=1,n
        x(i) = 1.d0/float(i) 
      end do   
      return                

      
12    continue
c                             Diagonal 2
      do i=1,n
        x(i) = 1.d0 
      end do   
      return     

      
13    continue
c                             Hager
      do i=1,n
        x(i) = 1.d0 
      end do   
      return     
      
14    continue
c                             Generalized Tridiagonal 1
      do i=1,n
        x(i) = 2.d0 
      end do   
      return     
      
15    continue
c                             Extended Tridiagonal 1
      do i=1,n
        x(i) = 2.d0 
      end do   
      return     
      
16    continue
c                             Extended Three Expo Terms
      do i=1,n
        x(i) = 0.1d0 
      end do   
      return      
      
17    continue                
c                             Generalized Tridiagonal 2
      do i=1,n
        x(i) = -1.d0 
      end do   
      return      
      
18    continue
c                             Diagonal 3 (1c1c)
      do i=1,n
        x(i) = 1.d0 
      end do   
      return     
      
19    continue
c                             Diagonal Full Borded
      do i=1,n
        x(i) = 0.1d0 
      end do   
      return    
      
20    continue
c                             Extended Himmelblau
      do i=1,n
        x(i) = 1.d0 
      end do   
      return 


21    continue
c                             Extended Powell
        i=1
921     x(i)  =  3.d0
        x(i+1)= -1.d0
        x(i+2)=  0.d0
        x(i+3)=  1.d0
        i=i+4
        if(i.le.n) go to 921             
      return


22	continue
c                             Tridiagonal Double Borded TR-DB1
      do i=1,n
        x(i) = -1.d0 
      end do   
      return 
   
     
23    continue
c                             Extended PSC1
        i=1
923     x(i)  =  3.d0
        x(i+1)=  0.1d0
        i=i+2
        if(i.le.n) go to 923       
      return
      
      
24    continue
c                             Extended BD1
      do i=1,n
        x(i) = 0.11d0 
      end do   
      return      
 
      
25    continue
c                             Extended Maratos
        i=1
925     x(i)  =  1.1d0
        x(i+1)=  0.1d0
        i=i+2
        if(i.le.n) go to 925       
      return        
 
      
26    continue
*                              FH1
*   Full Hessian FH1 (Summ of Quadratics, Quadratic inside)
*
      do i=1,n
        x(i) = float(i)/float(n) 
      end do   
      return


27    continue
c                             Extended Cliff
        i=1
927     x(i)  =   0.d0
        x(i+1)=  -0.1d0
        i=i+2
        if(i.le.n) go to 927       
      return 



28    continue
c                             Quadratic Diagonal Perturbed
      do i=1,n
        x(i) = 0.5d0 
      end do   
      return      
      

29	continue
*                              FH2
*    Full Hessian FH2 (Quadratic, perturbed with sin/cos)
*
      do i=1,n
        x(i) = 1.d0 
      end do
	return


30	continue
*                              FH3
*    Full Hessian FH3 (Quartic, perturbed with sin/cos)
*
      do i=1,n
        x(i) = 1.d0 
      end do
	return


31    continue
c                            TR-DBAD
*              Tridiagonal Double Borded Arrow-Down   
*                      
        i=1
931     x(i)  =  1.d0
        x(i+1)= -1.d0
        i=i+2
        if(i.le.n) go to 931      
      return 


32    continue
c                            TR-WHITEHOLST
*                   Tridiagonal White & Holst (c=4)
*                        
        i=1
932     x(i)  = -1.2d0
        x(i+1)=  1.0d0
        i=i+2
        if(i.le.n) go to 932 
	return


33    continue
c                            D-DBAUP3
*                 Diagonal Double Borded Arrow Up       
*                   
        i=1
933     x(i)  =  4.d0
        x(i+1)=  0.d0
        i=i+2
        if(i.le.n) go to 933      
      return 




34	continue
c                            TRIDIA
*                         Tridiagonal 
*
	do i=1,n
	  x(i) = 1.d0
	end do
	return
                             


35	continue
c                            ARWHEAD
*                Diagonal Bouble borded Arrow Down   
*
	do i=1,n
	  x(i) = 1.d0
	end do
	return


36    continue
c                            NONDIA (CUTE)
*                  Diagonal Bouble borded Arrow Up
      do i=1,n
        x(i) = -1.d0 
      end do   
      return  


37    continue
c                             Extended Wood   WOODS (CUTE)
        i=1
937     x(i)  =  -3.d0
        x(i+1)=  -1.d0
        x(i+2)=  -3.d0
        x(i+3)=  -1.d0
        i=i+4
        if(i.le.n) go to 937       
      return 


38    continue
c                             Extended Hiebert 
      do i=1,n
        x(i) = 0.00001d0 
      end do   
      return


39    continue
c                             BDQRTIC (CUTE)
      do i=1,n
        x(i) = 1.d0 
      end do   
      return 


40    continue
c                             DQDRTIC (CUTE)
      do i=1,n
        x(i) = 3.d0 
      end do   
      return  


41    continue
c                             EG2 (CUTE) 
      do i=1,n
        x(i) = 1.d0 
      end do   
      return  


42    continue
c                             EDENSCH (CUTE)
      do i=1,n
        x(i) = 0.d0 
      end do   
      return


43    continue
c                             Broyden Pentadiagonal (CUTE)
      do i=1,n
        x(i) = -1.d0 
      end do   
      return


44    continue
c                             Almost Perturbed Quadratic
      do i=1,n
        x(i) = 0.5d0 
      end do   
      return 


 
45    continue
c                             Almost Perturbed Quartic
      do i=1,n
        x(i) = 0.5d0 
      end do   
      return



46    continue
c                              FLETCHCR (CUTE)
      do i=1,n
        x(i) = 0.5d0 
      end do   
      return


47    continue
c                             ENGVAL1 (CUTE) 
      do i=1,n
        x(i) = 2.d0 
      end do   
      return



48    continue
c                             DENSCHNA (CUTE)
      do i=1,n
        x(i) = 8.d0 
      end do   
      return 



49    continue
c                             DENSCHNB (CUTE)
      do i=1,n
        x(i) = 1000.d0 
      end do   
      return



50    continue
c                             DENSCHNC  (CUTE)
      do i=1,n
        x(i) = 8.d0 
      end do   
      return



51    continue
c                             DENSCHNF (CUTE)      
        i=1
951     x(i)  =  100.d0
        x(i+1)= -100.d0
        i=i+2
        if(i.le.n) go to 951
      return



52    continue
c                             SINQUAD  (CUTE)
      do i=1,n
        x(i) = 0.1d0 
      end do   
      return  


53    continue
c                              HIMMELBG (CUTE)
      do i=1,n
        x(i) = 1.5d0 
      end do   
      return



54    continue
c                              HIMMELBH (CUTE)
      do i=1,n
        x(i) = 1.5d0 
      end do   
      return



55    continue
c                              DIXON3DQ (CUTE)
      do i=1,n
        x(i) = -1.d0 
      end do   
      return



56    continue
c                              BIGGSB1 (CUTE)
      do i=1,n
        x(i) = 0.d0            
      end do
      return



57    continue
c                              Perturbed Quadratic
      do i=1,n
        x(i) = 0.5d0 
      end do   
      return 


58    continue
c                              GENROSNB (CUTE)                       
        i=1
958     x(i)  = -1.2d0
        x(i+1)=  1.d0
        i=i+2
        if(i.le.n) go to 958
      return


59    continue      
c                              QP1 Extended Quadratic Penalty
      do i=1,n
        x(i) = 1.d0 
      end do   
      return 



60    continue
c                              QP2 Extended Quadratic Penalty
      do i=1,n
        x(i) = 1.d0 
      end do   
      return


61	continue
c                              STAIRCASE S1  
	do i=1,n
	  x(i) = 1.d0
	end do
	return




62	continue
c                              STAIRCASE S2  
	do i=1,n
	  x(i) = 1.d0
	end do
	return




63	continue
c                              S3  
	do i=1,n
	  x(i) = 2.d0
	end do
	return


64    continue
c                             Trigonometric        
      do i=1,n
        x(i) = 0.2d0 
      end do   
      return



65    continue
c                             QP3 Extended Quadratic Penalty       
      do i=1,n
        x(i) = 1.d0 
      end do   
      return



66    continue
c                             EG1       
      do i=1,n
        x(i) = 2.d0 
      end do   
      return



67    continue
c                             GENROSEN-2 
        i=1
967     x(i)  =  -1.2d0
        x(i+1)=   1.d0
        i=i+2
        if(i.le.n) go to 967      
      return                



68    continue
c                             PRODsin (m=n-1)
      do i=1,n
        x(i) = 5.d0 
      end do   
      return




69    continue
c                             PROD1 (m=n)      
      do i=1,n
        x(i) = 1.d0 
      end do   
      return





70    continue
c                             PRODcos (m=n-1)
      do i=1,n
        x(i) = 5.d0 
      end do   
      return




71    continue
c                             PROD2 (m=1)      
      do i=1,n
        x(i) = 15.d0 
      end do   
      return



72    continue
c                             ARGLINB (m=5) 
        i=1
972     x(i)  =  0.01d0
        x(i+1)=  0.001d0
        i=i+2
        if(i.le.n) go to 972      
      return  



73    continue
c                             DIXMAANA (CUTE)
      do i=1,n
        x(i) = 2.d0 
      end do   
      return 



74    continue
c                             DIXMAANB (CUTE)
      do i=1,n
        x(i) = 2.d0 
      end do   
      return



75    continue
c                             DIXMAANC (CUTE)
      do i=1,n
        x(i) = 2.d0 
      end do   
      return



76    continue
c                             DIXMAAND (CUTE)
      do i=1,n
        x(i) = 2.d0 
      end do   
      return




77    continue
c                             DIXMAANE (CUTE)
      do i=1,n
        x(i) = 2.d0 
      end do   
      return



78    continue
c                             VARDIM (CUTE)
      do i=1,n
        x(i) = 1.d0-float(i)/float(n) 
      end do   
      return



79    continue
c                             DIAG-AUP1
      do i=1,n
        x(i) = 4.d0 
      end do   
      return


80    continue
c                             ENGVAL8
      do i=1,n
        x(i) = 2.d0 
      end do   
      return




***
      end  
            

c------------------------------------------------ End INIPOINT













***************************************************************
*                            Date created: October 28, 2004   *
*                            80 problems                      *
*                                                             *
*                                                             *
*       TEST FUNCTIONS FOR UNCONSTRAINED OPTIMIZATION         *
*      ===============================================        *
*                                                             *
* Subroutine for functions and their gradient specification.  *
* This is a user subroutine:                                  *
*                                                             *
*  The calling sequence is:                                   *
*                                                             *
*   call evalfg(n,x,f,g, nexp)                                *
* where:                                                      *
*   n    (integer) the number of variables,                   *
*   x    (double)  array with the initial point,              *
*   f    (double)  function value in point x,                 *
*   g    (double)  array with gradient value in point x,      *
*   nexp (integer) parameter specifying the number of the     *
*                  problem considered in a train of           *
*                  experiments.                               *
*                                                             *  
* ----------------------------------------------------------- *
*                                                             *
* 57 problems:  October 28, 2004                              *
* 66 problems:  March 29,   2005                              *
* 70 problems:  April 28,   2005                              *
* 75 problems:  April 19,   2006                              *
* 80 problems:  January 28, 2013                              *
*                                                             *  
*                                               Neculai Andrei*
***************************************************************    
*                                                             
      subroutine evalfg(n,x,f,g, nexp)
      real*8 x(n), f, g(n)  
      
      real*8 t1,t2,t3,t4, c, d
      real*8 s, temp(1000000), temp1, tsum, sum
      real*8 u(1000000), v(1000000), t(1000000) 
      real*8 u1, v1, c1, c2
      real*8 alpha, beta, gamma, delta  
                       
      integer k1, k2, k3, k4      
                         
*      
      go to ( 1, 2, 3, 4, 5, 6, 7, 8, 9,10,
     *       11,12,13,14,15,16,17,18,19,20,
     *       21,22,23,24,25,26,27,28,29,30,
     *       31,32,33,34,35,36,37,38,39,40,
     *       41,42,43,44,45,46,47,48,49,50,
     *       51,52,53,54,55,56,57,58,59,60,
     *       61,62,63,64,65,66,67,68,69,70,
     *       71,72,73,74,75,76,77,78,79,80) nexp  
     
     

   
cF1                                FREUROTH (CUTE)
*                          Extended Freudenstein & Roth 
*
*                      Initial Point: [0.5, -2, ...,0.5, -2].
*

1     continue
      f = 0.d0
      j=1         
      
      do i=1,n/2 
        t1=-13.d0+x(2*i-1)+5.d0*x(2*i)*x(2*i)-x(2*i)**3-2.d0*x(2*i) 
        t2=-29.d0+x(2*i-1)+x(2*i)**3+x(2*i)**2-14.d0*x(2*i)   
        
        f = f + t1*t1 + t2*t2
*--        
        g(j)  =2.d0*(t1+t2)   
        g(j+1)=2.d0*t1*(10.d0*x(2*i)-3.d0*x(2*i)*x(2*i)-2.d0) +
     *         2.d0*t2*(3.d0*x(2*i)*x(2*i)+2.d0*x(2*i)-14.d0)   
        j=j+2        
      end do   
      
      return

             


cF2                           Extended Trigonometric ET1
*
*                          Initial Point: [0.2, 0.2, ....,0.2].

2     continue      
      s= float(n)
      do i=1,n
        s = s - dcos(x(i))
      end do  
	 
      do i=1,n
        temp(i) = s + float(i)*(1.d0-dcos(x(i))) - dsin(x(i))
      end do 

      f = 0.d0
      do i=1,n
        f = f + temp(i)**2
      end do   
*--
      s=0.d0
      do i=1,n
        s = s + temp(i)
      end do
      
      do i=1,n
        g(i) = 2.d0*s*dsin(x(i)) + 
     +            2.d0*temp(i)*(float(i)*dsin(x(i))-dcos(x(i)))  
      end do      
      return





 
cF3                          SROSENBR (CUTE)
*                     Extended Rosenbrock function   
*
* Initial point: [-1.2, 1, -1.2, 1, ..........., -1.2, 1]      
      
3     continue                           
      c=1000.d0                           
                           
      f=0.d0
      do i=1,n/2
        f = f + c*(x(2*i)-x(2*i-1)**2)**2 + (1.d0-x(2*i-1))**2
      end do  
*--               
      j=1
      do i=1,n/2
        g(j)   = -4.d0*c*x(2*i-1)*(x(2*i)-x(2*i-1)**2) - 
     *            2.d0*(1.d0-x(2*i-1))
        g(j+1) =  2.d0*c*(x(2*i)-x(2*i-1)**2)
        j = j + 2
      end do
      
      return








cF4                    Extended White & Holst function
*
*           Initial point: [-1.2, 1, -1.2, 1, ..........., -1.2, 1]      
      
4     continue                           
      c=1.d0                           
                           
      f=0.d0
      do i=1,n/2
        f = f + c*(x(2*i)-x(2*i-1)**3)**2 + (1.d0-x(2*i-1))**2
      end do  
*--               
      j=1
      do i=1,n/2
        g(j)   = -6.d0*c*x(2*i-1)*x(2*i-1)*(x(2*i)-x(2*i-1)**3) - 
     *            2.d0*(1.d0-x(2*i-1))
        g(j+1) =  2.d0*c*(x(2*i)-x(2*i-1)**3)
        j = j + 2
      end do
      
      return






*
cF5                  Extended Beale Function  BEALE (CUTE)
*
*                   Initial Point: [1, 0.8, ...., 1, 0.8]
*
5     continue
*      
      f=0.d0

      j=1      
      do i=1,n/2 
        t1=1.5d0  -x(2*i-1)+x(2*i-1)*x(2*i)
        t2=2.25d0 -x(2*i-1)+x(2*i-1)*x(2*i)*x(2*i)
        t3=2.625d0-x(2*i-1)+x(2*i-1)*x(2*i)*x(2*i)*x(2*i)   
        
        f = f + t1*t1 + t2*t2 + t3*t3
*--
        g(j)  =2.d0*t1*(-1.d0+x(2*i)) +
     *         2.d0*t2*(-1.d0+x(2*i)**2) +
     *         2.d0*t3*(-1.d0+x(2*i)**3)   
        g(j+1)=2.d0*t1*x(2*i-1) +
     *         2.d0*t2*2.d0*x(2*i-1)*x(2*i) +
     *         2.d0*t3*3.d0*x(2*i-1)*x(2*i)*x(2*i)   
        j=j+2        
      end do   
*   
      return



cF6                   Extended Penalty Function  U52 (MatrixRom)
*
*                           Intial Point: [1,2,3,.....,n].
*
6     continue      
      temp1=0.d0
      do i=1,n
        temp1 = temp1 + x(i)**2
      end do
*      
      f = (temp1 - 0.25d0)**2
      
      do i=1,n-1
        f = f + (x(i)-1.d0)**2 
      end do
*--        
      do i=1,n-1
        g(i) = 2.d0*(x(i)-1.d0) + 4.d0*x(i)*(temp1-0.25d0)
      end do                                              
      
      g(n) = 4.d0*x(n)*(temp1-0.25d0)
      return


cF7                        Perturbed Quadratic function
*
*                        Initial Point:  [1, 1, ......,1].
*     
7     continue
      
      temp1 = 0.d0
      do i=1,n
        temp1 = temp1 + x(i)
      end do  
      
      f = temp1*temp1/100.d0
      
      do i=1,n
        f = f + float(i)*x(i)**2      
      end do
*--        
      do i=1,n
        g(i) = float(i) * 2.d0 * x(i) + temp1/50.d0
      end do        
      return



cF8                        Raydan 1 Function
*
*                   Initial point: [1, 1, ..., 1]

8     continue
      
      f=0.d0
      do i=1,n
        f = f + float(i) * (dexp(x(i))-x(i)) / 10.d0
      end do  
*--        
      do i=1,n
        g(i) = float(i) * (dexp(x(i)) - 1.d0) / 10.d0      
      end do  
      return




cF9                             Raydan 2 Function
*
*                       Initial Point: [1, 1, .....,1]
*  
9     continue
      
      f=0.d0
      do i=1,n
        f = f + dexp(x(i)) - x(i)
      end do
*--      
      do i=1,n
        g(i) = dexp(x(i)) - 1.d0
      end do
      return




cF10                             TR-Summ of quadratics Function
*
*                               Initial Point: [1, 1, ....., 1]
*  
10    continue
      
	c = 100000.d0
      f = 0.d0  

      do i=1,n-1
        f = f + x(i)*x(i) + c*(x(i+1)+x(i)*x(i))**2
      end do
*--
      g(1) = 2.d0*x(1) + 4.d0*c*(x(2)+x(1)*x(1))*x(1)
	        
      do i=2,n-1
        g(i) = 2.d0*c*(x(i)+x(i-1)**2) + 2.d0*x(i) +
     *         4.d0*c*(x(i+1)+x(i)*x(i))*x(i)   
      end do

      g(n) = 2.d0*c*(x(n)+x(n-1)**2)
      
      return



cF11                             Diagonal1 Function
*
*                       Initial Point: [1/1, 1/2, ....., 1/n]
*  
11    continue
      
      f=0.d0
      
      do i=1,n
        f = f + dexp(x(i)) - x(i)/float(i)
      end do
*--      
      do i=1,n
        g(i) = dexp(x(i)) - 1.d0/float(i)
      end do
      return



cF12                         Diagonal2 Function
*
*                        Initial Point: [1,1,...,1] 

12    continue
      
      f=0.d0
      do i=1,n
        f = f + dexp(x(i)) - float(i)*dsin(x(i))
      end do
*--      
      do i=1,n
        g(i) = dexp(x(i)) - float(i)*dcos(x(i))
      end do
      return



cF13                          Hager Function
*
*                        Initial Point: [1,1,...,1] 


13    continue
      f=0.d0
      do i=1,n
        f = f + dexp(x(i)) - x(i)*sqrt(float(i))
      end do
*--      
      do i=1,n
        g(i) = dexp(x(i)) - sqrt(float(i))
      end do
      return



cF14                 Generalized Tridiagonal-1 Function
*
*                        Initial Point: [2,2,...,2] 


14    continue
      
      f=0.d0
      do i=1,n-1
        u(i) = x(i) + x(i+1) - 3.d0
        v(i) = x(i) - x(i+1) + 1.d0
      end do

      do i=1,n-1
        f = f + u(i)**2 + v(i)**4
      end do
*--                       
      g(1) = 2.d0*u(1) + 4.d0*v(1)**3

      do i=2,n-1
        g(i) = 2.d0*u(i-1) - 4.d0*v(i-1)**3 + 2.d0*u(i) + 4.d0*v(i)**3
      end do

      g(n) = 2.d0*u(n-1) - 4.d0*v(n-1)**3                         

      return



cF15                   Extended Tridiagonal-1 Function
*
*                        Initial Point: [2,2,...,2] 


15    continue
      
      f=0.d0
      do i=1,n/2
        u(i) = x(2*i-1) + x(2*i) - 3.d0
        v(i) = x(2*i-1) - x(2*i) + 1.d0
      end do
      do i=1,n/2
        f = f + u(i)**2 + v(i)**4
      end do
*--                       
      j=1
      do i=1,n/2
        g(j)   = 2.d0*u(i) + 4.d0*v(i)**3
        g(j+1) = 2.d0*u(i) - 4.d0*v(i)**3
        j=j+2
      end do                        

      return



cF16                        Extended Three Exponential Terms
*
*                          Intial Point: [0.1,0.1,.....,0.1].
*

16    continue
      f=0.d0   

      j=1
      do i=1,n/2
        t1= x(2*i-1) + 3.d0*x(2*i) - 0.1d0
        t2= x(2*i-1) - 3.d0*x(2*i) - 0.1d0
        t3=-x(2*i-1) - 0.1d0 

        f = f + dexp(t1) + dexp(t2) + dexp(t3)    
*--        
        g(j)   = dexp(t1) + dexp(t2) - dexp(t3)
        g(j+1) = 3.d0*dexp(t1) - 3.d0*dexp(t2)        
        j=j+2
      end do  
      return




cF17                       Generalized Tridiagonal-2  
*                                Penta Diagonal
*                    Initial point: [-1, -1, .........., -1., -1] 
 
17    continue
      f = 0.d0
      u(1) = 5.d0*x(1)-3.d0*x(1)**2-x(1)**3-3.d0*x(2)+1.d0
      
      do i=2,n-1
        u(i)=5.d0*x(i)-3.d0*x(i)**2-x(i)**3-x(i-1)-3.d0*x(i+1)+1.d0
      end do      
                
      u(n)=5.d0*x(n)-3.d0*x(n)**2-x(n)**3-x(n-1)+1.d0        
      
      do i=1,n   
        f = f + u(i)**2
        v(i) = 5.d0 -6.d0*x(i) -3.d0*x(i)**2
      end do    
      
*--

      g(1) = 2.d0*u(1)*v(1) - 2.d0*u(2)
      
      do i=2,n-1
        g(i) = -6.d0*u(i-1) + 2.d0*u(i)*v(i) - 2.d0*u(i+1)
      end do
      
      g(n) = -6.d0*u(n-1) + 2.d0*u(n)*v(n)  
      return



cF18                       Diagonal3 (1c1c) Function  
*                      
*                Initial point: [1, 1, .........., 1., 1] 
 
18    continue
      
      c=10000.d0

      f = 0.d0
      do i=1,n/2
        f = f + (x(2*i-1)**2 + c*x(2*i)**2)/2.d0
      end do                  
*--
      j=1
      do i=1,n/2
        g(j)   =   x(2*i-1)
        g(j+1) = c*x(2*i)
        j=j+2
      end do  
      return


cF19                       Diagonal Full Borded
*                      
*                Initial point: [0.1, 0.1, .........., 0.1] 
 
19    continue
      f=(x(1)-1.d0)**4 + (x(n)**2-x(1)**2)**2
      
      do i=1,n-2
        temp(i) = sin(x(i+1)-x(n)) - x(1)**2 - x(i+1)**2
        f = f + temp(i)*temp(i)
      end do
*--      
      g(1) = 4.d0*(x(1)-1.d0)**3 - 4.d0*x(1)*(x(n)**2-x(1)**2) 
      do i=1,n-2
        g(1) = g(1) - 4.d0*temp(i)*x(1)
      end do
      
      do i=2,n-1
        g(i) = 2.d0*temp(i-1)*(cos(x(i)-x(n))-2.d0*x(i))
      end do
      
      g(n) = 4.d0*x(n)*(x(n)**2-x(1)**2)
      do i=1,n-2
        g(n) = g(n) - 2.d0*temp(i)*cos(x(i+1)-x(n))    
      end do
      return



cF20                             HIMMELBC (CUTE)
*                         Extended Himmelblau Function
*
*                           Initial Point: [1, 1, ....., 1]

20    continue
      
      f=0.d0

      j=1
      do i=1,n/2
        u1 = x(2*i-1)**2 + x(2*i)    - 11.d0
        v1 = x(2*i-1)    + x(2*i)**2 - 7.d0 
        
        f = f + u1*u1 + v1*v1
*--
        g(j)   = 4.d0*u1*x(2*i-1) + 2.d0*v1
        g(j+1) = 2.d0*u1          + 4.d0*v1*x(2*i)
        j=j+2
      end do

      return



cF21                              Extended Powell 
*
*                      Initial Point: [3, -1, 0, 1, ......].
*-------------------------------------------------------------------
*
21    continue
      
      f=0.d0
      
      j=1    
      do i=1,n/4  
        t1= x(4*i-3) + 10.d0*x(4*i-2) 
        t2= x(4*i-1) - x(4*i)
        t3= x(4*i-2) - 2.d0*x(4*i-1)
        t4= x(4*i-3) - x(4*i)     
       
        f = f + t1*t1 + 5.d0*t2*t2 + t3**4 + 10.d0*t4**4      
*--        
        g(j)  =   2.d0*t1 + 40.d0*t4**3
        g(j+1)=  20.d0*t1 +  4.d0*t3**3
        g(j+2)=  10.d0*t2 -  8.d0*t3**3
        g(j+3)= -10.d0*t2 - 40.d0*t4**3
        
        j=j+4        
      end do               
      return




cF22                        TR-DB1                                 
*                   Tridiagonal Double Borded 
*               Initial Point:  [-1, -1, ...., -1]
*
22    continue
	f = (x(1)-1.d0)**2

	do i=1,n-1
	  temp(i) = x(1) - 0.5d0*x(i)**2 - 0.5d0*x(i+1)**2
	  f = f + temp(i)*temp(i)
      end do
*--
	g(1) = 2.d0*(x(1)-1.d0) + 2.d0*temp(1)*(1.d0-x(1))

	do i=2,n-1
	  g(1) = g(1) + 2.d0*temp(i)
	end do
	
	do i=2,n-1
	  g(i) = -2.d0*x(i)*(temp(i-1) + temp(i))
	end do

	g(n) = -2.d0*x(n)*temp(n-1)    	
	
      return






cF23                       Extended PSC1 Function
*
*                   Initial point: [3, 0.1, ..., 3, 0.1]

23    continue 

      f = 0.d0
      do i=1,n/2
        f = f + (x(2*i-1)**2 +x(2*i)**2 +x(2*i-1)*x(2*i))**2 
     *        + (dsin(x(2*i-1)))**2 + (dcos(x(2*i)))**2    
      end do
*--     
      j=1
      do i=1,n/2
      g(j)   = 2.d0*(x(2*i-1)**2+x(2*i)**2+x(2*i-1)*x(2*i)) *
     +              (2.d0*x(2*i-1)+x(2*i)) +
     +              2.d0*(dsin(x(2*i-1)))*(dcos(x(2*i-1)))
     
      g(j+1) = 2.d0*(x(2*i-1)**2+x(2*i)**2+x(2*i-1)*x(2*i)) *
     +              (2.d0*x(2*i)+x(2*i-1)) -
     +              2.d0*(dcos(x(2*i)))*(dsin(x(2*i)))
      j=j+2
      end do
      return






cF24                  Extended Block Diagonal BD1 Function 
*
*                     Initial Point: [0.1, 0.1, ..., 0.1].
*
24    continue
      
      f = 0.d0
      
      j=1
      do i=1,n/2
        t1 = x(2*i-1)**2 + x(2*i)**2 - 2.d0
        t2 = dexp(x(2*i-1)) - x(2*i)
        
        f = f + t1*t1 + t2*t2
*--        
        g(j)   = 4.d0*t1*x(2*i-1) + 2.d0*t2*dexp(x(2*i-1))
        g(j+1) = 4.d0*t1*x(2*i) - 2.d0*t2
        j=j+2
      end do
      return



cF25                           Extended Maratos Function 
*
*                      Initial Point: [1.1, 0.1, ...,1.1, 0.1].
*
25    continue
      
      c = 1.d0
      f = 0.d0
      
      j=1
      do i=1,n/2
        t1 = x(2*i-1)**2 + x(2*i)**2 - 1.d0
        
        f = f + (x(2*i-1) + c*t1*t1)
*--        
        g(j)   = 1.d0 + 4.d0 * c * t1 * x(2*i-1)
        g(j+1) =        4.d0 * c * t1 * x(2*i)
        j=j+2
      end do
      return



cF26                         FH (m=500)                               
*        Full Hessian FH1 (Summ of Quadratics, Quadratic inside)
*
26    continue
      m=500
      f=0.d0         
      
      do i=1,m
        u(i)=0.d0
        do j=1,n
          u(i) = u(i) + float(i)*float(j)*x(j)*x(j)
        end do  

        f = f + (u(i)-1.d0)**2
      end do  
*--
      do j=1,n
        g(j) = 0.d0
        do i=1,m
          g(j) = g(j) + 4.d0*(u(i)-1.d0)*float(i)*float(j)*x(j)
        end do
      end do                              

	return


cF27                             Extended CLIFF (CUTE)                
*
*                       Initial Point:  [0, -0.1, ......, 0, -0.1].
*                     

27    continue
      
      f=0.d0
        
      j=1  
      do i=1,n/2
        temp1 = (x(2*i-1)-3.d0)/100.d0
        
      f = f+temp1*temp1-(x(2*i-1)-x(2*i))+dexp(20.d0*(x(2*i-1)-x(2*i)))
*--        
       g(j)   = temp1/50.d0 - 1.d0 + 20.d0*dexp(20.d0*(x(2*i-1)-x(2*i)))
       g(j+1) = 1.d0 - 20.d0*dexp(20.d0*(x(2*i-1)-x(2*i)))
       j=j+2
      end do
      
      return   





cF28                   Quadratic Diagonal Perturbed Function
*
*                     Initial Point:  [0.5, 0.5, ......, 0.5].
*     
28    continue
      
      temp1 = 0.d0
      do i=1,n
        temp1 = temp1 + x(i)
      end do  
      
      f = temp1*temp1
      
      do i=1,n
        f = f + (float(i)/100.d0) * x(i)**2      
      end do
*--     
      do i=1,n
        g(i) = float(i) * x(i) / 50.d0 + 2.d0*temp1
      end do        
      return




cF29                        FH2                                
*      Full Hessian FH2 (Quadratic, perturbed with sin/cos)
*           Initial Point:  [1, 1, ......, 1].
*
29    continue
	s = 0.d0
	do i=1,n
	  s = s + x(i)
      end do
	
	f = s*s

	do i=1,n
	  f = f + float(i)*(dsin(x(i)) + dcos(x(i)))/1000.d0
      end do	
*--
	do i=1,n
	  g(i) = 2.d0*s + float(i)*(dcos(x(i)) - dsin(x(i)))/1000.d0
	end do

	return




cF30                        FH3                                
*      Full Hessian FH3 (Quartic, perturbed with sin/cos)
*           Initial Point:  [1, 1, ......, 1].
*
30    continue
	s = 0.d0
	do i=1,n
	  s = s + x(i)**2
      end do
	
	f = s*s

	do i=1,n
	  f = f + float(i)*(dsin(x(i)) + dcos(x(i)))/1000.d0
      end do	
*--
	do i=1,n
	  g(i) = 4.d0*s*x(i) + float(i)*(dcos(x(i)) - dsin(x(i)))/1000.d0
	end do

	return



cF31                         NONDQUAR                    
*              Tridiagonal Double Borded Arrow-Down
*
31    continue

      f = (x(1)-x(2))**2 + (x(n-1)+x(n))**2      
      
      do i=1,n-2
        f = f + (x(i)+x(i+1)+x(n))**4
      end do
*--
      g(1) = 2.d0*(x(1)-x(2))+4.d0*(x(1)+x(2)+x(n))**3
      g(2) =-2.d0*(x(1)-x(2))+4.d0*(x(1)+x(2)+x(n))**3 +
     *                        4.d0*(x(2)+x(3)+x(n))**3
      
      do i=3,n-2
        g(i) = 4.d0*(x(i-1)+x(i)+x(n))**3 +
     *         4.d0*(x(i)+x(i+1)+x(n))**3
      end do          
      
      g(n-1) = 4.d0*(x(n-2)+x(n-1)+x(n))**3 +
     *         2.d0*(x(n-1)+x(n))
      
      g(n) = 2.d0*(x(n-1)+x(n))
      do i=1,n-2
        g(n) = g(n) + 4.d0*(x(i)+x(i+1)+x(n))**3
      end do  
              
      return



cF32                         TR-WHITEHOLST
*                    Tridiagonal. White-Holst (c=4)
*             Initial point x =[-1.2,  1, ..., -1.2, 1]
*
32    continue
      c = 4.d0
      f = 0.d0
      
      do i=1,n-1
        f = f + c*(x(i+1)-x(i)**3)**2 + (1.d0-x(i))**2
      end do
*--
      g(1) = -6.d0*c*(x(2)-x(1)**3)*x(1)*x(1) - 2.d0*(1.d0-x(1))
      
      do i=2,n-1
        g(i) = 2.d0*c*(x(i)-x(i-1)**3) -
     *         6.d0*c*(x(i+1)-x(i)**3)*x(i)*x(i) -
     *         2.d0*(1.d0-x(i))
      end do
      
      g(n) = 2.d0*c*(x(n)-x(n-1)**3)
      
      return



cF33                         D-DBAUP3                              
*                 Diagonal Double Borded Arrow Up
*               Initial point: x0=[4, 0, ....,4,0]
*
33    continue
      f=0.d0

      do i=1,n
        f = f + 4.d0*(x(i)*x(i) - x(1))**2 + (x(i)-1.d0)**2
      end do  
*--
      g(1) = 2.d0*(x(1)-1.d0) + 8.d0*(x(1)*x(1)-x(1))*(2.d0*x(1)-1.d0)
      do i=2,n
        g(1) = g(1) - 8.d0*(x(i)*x(i)-x(1))
      end do
      
      do i=2,n
        g(i) = 16.d0*x(i)*(x(i)*x(i)-x(1)) + 2.d0*(x(i)-1.d0)
      end do  
        
      return




cF34                           TRIDIA  (CUTE)                             
*
*                     Initial point x0=[1,1,...,1].
*
*
34    continue
*                
      alpha=5.d0
      beta =1.d0
      gamma=1.d0
      delta=1.d0
        
      f=gamma*(delta*x(1)-1.d0)**2
      
      do i=2,n
        f = f + float(i)*(alpha*x(i)-beta*x(i-1))**2
      end do
*--
      g(1) = 2.d0*gamma*(delta*x(1)-1.d0)*delta - 
     *       4.d0*(alpha*x(2)-beta*x(1))*beta
     
      do i=2,n-1
        g(i) = 2.d0*float(i)*(alpha*x(i)-beta*x(i-1))*alpha -
     *         2.d0*float(i+1)*(alpha*x(i+1)-beta*x(i))*beta
      end do         

      g(n) = 2.d0*float(n)*(alpha*x(n)-beta*x(n-1))*alpha

      return




cF35                       ARWHEAD  (CUTE)                             
*                          
*                     Initial point x0=[1,1,...,1].
*
*
35    continue

      f=0.d0
      do i=1,n-1
        f = f + (-4.d0*x(i)+3.d0) + (x(i)**2+x(n)**2)**2
      end do
*--      
      do i=1,n-1
        g(i) = -4.d0 + 4.d0*x(i)*(x(i)**2+x(n)**2)
      end do
      
      g(n) = 0.d0
      do i=1,n-1
        g(n) = g(n) + 4.d0*x(n)*(x(i)**2+x(n)**2)
      end do                    
      
      return




cF36                                                         
*                       NONDIA (Shanno-78)  (CUTE)
*
*                     Initial point x0=[-1,-1,...,-1].
*
*
36    continue

      c=100.d0
      
      f=(x(1)-1.d0)**2 + c*(x(1)-x(1)**2)**2
      
      do i=2,n
        f = f + c*(x(1)-x(i)**2)**2
      end do
*--
      g(1)=2.d0*(x(1)-1.d0) + 2.d0*c*(x(1)-x(1)**2)*(1.d0-2.d0*x(1))        
      do i=2,n
        g(1) = g(1) + 2.d0*c*(x(1)-x(i)**2)
      end do  
      
      do i=2,n
        g(i) = -4.d0*c*x(i)*(x(1)-x(i)**2)
      end do        
              
      return



cF37                           Extended Wood Function
*                                   WOODS (CUTE)
*
*                         Initial Point: [-3,-1,-3,-1,......]
*
37    continue
      
      f=0.d0
      
      j=1
      do i=1,n/4
        f = f + 100.d0*(x(4*i-3)**2-x(4*i-2))**2 
     *        +        (x(4*i-3)-1.d0)**2  
     *        +  90.d0*(x(4*i-1)**2-x(4*i))**2            
     *        +        (1.d0-x(4*i-1))**2
     *        +  10.1d0*(x(4*i-2)-1.d0)**2
     *        +  10.1d0*(x(4*i)  -1.d0)**2
     *        +  19.8d0*(x(4*i-2)-1.d0)*(x(4*i)-1.d0)
*--     
        g(j)   = 400.d0*(x(4*i-3)**2-x(4*i-2))*x(4*i-3)
     *             + 2.d0*(x(4*i-3)-1.d0)
        g(j+1) =-200.d0*(x(4*i-3)**2-x(4*i-2))
     *             + 20.2d0*(x(4*i-2)-1.d0)   
     *             + 19.8d0*(x(4*i)-1.d0)
        g(j+2) = 360.d0*(x(4*i-1)**2-x(4*i))*x(4*i-1)
     *             - 2.d0*(1.d0-x(4*i-1))   
        g(j+3) =-180.d0*(x(4*i-1)**2-x(4*i))
     *             + 20.2d0*(x(4*i)  -1.d0)
     *             + 19.8d0*(x(4*i-2)-1.d0)   
        j=j+4
      end do
      return



cF38                          Extended Hiebert Function
*
*                             Initial Point: [0,0,...0].

38    continue
      
      c1 = 10.d0
      c2 = 500.d0
      
      f = 0.d0
      
      j=1
      do i=1,n/2
        f = f + (x(2*i-1)-c1)**2 + (x(2*i-1)*x(2*i)-c2)**2
*--
        g(j)   =   2.d0*(x(2*i-1)-c1) 
     *           + 2.d0*(x(2*i-1)*x(2*i)-c2)*x(2*i)       
        g(j+1) =   2.d0*(x(2*i-1)*x(2*i)-c2)*x(2*i-1)
        j=j+2
      end do
      return



cF39                           BDQRTIC (CUTE)
*
*                     Initial point x0=[1.,1.,...,1.].
*
*
39    continue
*        
      n4=n-4
      f=0.d0
        
      do i=1,n4
        temp(i) = x(i)**2 + 2.d0*x(i+1)**2 + 3.d0*x(i+2)**2 
     *                    + 4.d0*x(i+3)**2 + 5.d0*x(n)**2
      end do           
        
      do i=1,n4
        f = f + (-4.d0*x(i)+3.d0)**2 + temp(i)**2
      end do                                  
*--
      g(1) = -8.d0*(-4.d0*x(1)+3.d0) + 
     *        (4.d0*temp(1))*x(1)        
      g(2) = -8.d0*(-4.d0*x(2)+3.d0) + 
     *        (8.d0*temp(1)+  4.d0*temp(2))*x(2)
      g(3) = -8.d0*(-4.d0*x(3)+3.d0) + 
     *        (12.d0*temp(1)+ 8.d0*temp(2)+ 4.d0*temp(3))*x(3)
      g(4) = -8.d0*(-4.d0*x(4)+3.d0) + 
     *        (16.d0*temp(1)+12.d0*temp(2) +8.d0*temp(3)+
     *          4.d0*temp(4))*x(4)
        
      do i=5,n4
        g(i) = -8.d0*(-4.d0*x(i)+3.d0) + 
     *         (16.d0*temp(i-3)+12.d0*temp(i-2)+
     *           8.d0*temp(i-1)+4.d0*temp(i))*x(i)                                                        
      end do
        
      g(n4+1) =(16.d0*temp(n4-2)+12.d0*temp(n4-1)+8.d0*temp(n4))*x(n4+1)  
      g(n4+2) =(16.d0*temp(n4-1)+12.d0*temp(n4))*x(n4+2)
      g(n4+3) =(16.d0*temp(n4))*x(n4+3)
          
      tsum=0.d0
      do i=1,n4
        tsum = tsum + temp(i)
      end do
      g(n) = 20.d0*tsum*x(n)
             
      return




cF40                          DQDRTIC  (CUTE)                            
*                             
*                     Initial point x0=[3,3,3...,3].
*
*
40    continue
      
      c=1000.d0
      d=1000.d0
      
      f=0.d0
      do i=1,n-2
        f = f + (x(i)**2 + c*x(i+1)**2 + d*x(i+2)**2)
      end do
*--
      g(1) = 2.d0*x(1)
      g(2) = 2.d0*c*x(2) + 2.d0*x(2)
      
      do i=3,n-2
        g(i) = 2.d0*(1.d0+d+c)*x(i)
      end do
      
      g(n-1) = 2.d0*(c+d)*x(n-1)          
      g(n) = 2.d0*d*x(n)
      
      return




cF41                          EG2 (CUTE)                             
*                             
*                     Initial point x0=[1,1,1...,1].
*
*
41    continue

      f=0.5d0*dsin(x(n)*x(n))
      do i=1,n-1
        f = f + dsin(x(1)+x(i)*x(i)-1.d0)
      end do
*--
      g(1)=(1.d0+2.d0*x(1))*dcos(x(1)+x(1)*x(1)-1.d0)
      do i=2,n-1
        g(1) = g(1) + dcos(x(1)+x(i)*x(i)-1.d0)
      end do               
      
      do i=2,n-1
        g(i) = 2.d0*x(i)*dcos(x(1)+x(i)*x(i)-1.d0)
      end do
      
      g(n) = x(n)*dcos(x(n)*x(n))

      return




cF42                       EDENSCH Function  (CUTE)
*
*                      Initial Point: [0., 0., ..., 0.].


42    continue      

      f = 16.d0

      do i=1,n-1
        f = f + (x(i)-2.d0)**4 + 
     *          (x(i)*x(i+1)-2.d0*x(i+1))**2 +
     *          (x(i+1)+1.d0)**2
      end do
          
*--          
      g(1) = 4.d0*(x(1)-2.d0)**3 + 2.d0*x(2)*(x(1)*x(2)-2.d0*x(2))
      
      do i=2,n-1               
        g(i) = 2.d0*(x(i-1)*x(i)-2.d0*x(i))*(x(i-1)-2.d0) +
     *         2.d0*(x(i)+1.d0) +
     *         4.d0*(x(i)-2.d0)**3 +
     *         2.d0*x(i+1)*(x(i)*x(i+1)-2.d0*x(i+1))   
      end do                                        
      
      g(n) = 2.d0*(x(n-1)*x(n)-2.d0*x(n))*(x(n-1)-2.d0) +
     *       2.d0*(x(n)+1.d0) 
      
      return




cF43                       Broyden Pentadiagonal                              
*                             
*                  Initial point x0=[-1., -1., ..., -1.].
*
43    continue
*

      temp(1) = 3.d0*x(1) - 2.d0*x(1)*x(1)
      do i=2,n-1
        temp(i) = 3.d0*x(i)-2.d0*x(i)*x(i)-x(i-1)-2.d0*x(i+1)+1.d0
      end do
        
      temp(n) = 3.d0*x(n)-2.d0*x(n)*x(n)-x(n-1)+1.d0
      
      f = 0.d0
      
      do i=1,n
        f = f + temp(i)*temp(i)
      end do
*--
      g(1) = 2.d0*temp(1)*(3.d0-4.d0*x(1)) - 2.d0*temp(2)
      
      g(2) = 2.d0*temp(2)*(3.d0-4.d0*x(2)) - 2.d0*temp(3)
      
      do i=3,n-1
        g(i) = -4.d0*temp(i-1) 
     *         +2.d0*temp(i)*(3.d0-4.d0*x(i)) 
     *         -2.d0*temp(i+1)             
      end do
      
      g(n) = -4.d0*temp(n-1) + 2.d0*temp(n)*(3.d0-4.d0*x(n))
      
      return






cF44                   Almost Perturbed Quadratic                              
*                             
*                  Initial point x0=[0.5, 0.5, ...,0.5].
*
44    continue
*    

      f = ((x(1)+x(n))**2)/100.d0
      
      do i=1,n
        f = f + float(i)*x(i)*x(i)
      end do     
*--
      g(1) = 2.d0*x(1) + (x(1)+x(n))/50.d0
      
      do i=2,n-1
        g(i) = 2.d0*float(i)*x(i)
      end do
      
      g(n) = 2.d0*float(n)*x(n) + (x(1)+x(n))/50.d0
      
      return




cF45                   Almost Perturbed Quartic                              
*                             
*                  Initial point x0=[0.5, 0.5, ...,0.5].
*
45    continue
*    

      f = ((x(1)+x(n))**2)/100.d0
      
      do i=1,n
        f = f + float(i)*x(i)**4
      end do     
*--
      g(1) = 4.d0*x(1)**3 + (x(1)+x(n))/50.d0
      
      do i=2,n-1
        g(i) = 4.d0*float(i)*x(i)**3
      end do
      
      g(n) = 4.d0*float(n)*x(n)**3 + (x(1)+x(n))/50.d0
      
      return





cF46                       FLETCHCR (CUTE)              
*
*                  Initial Point: [0.5,0.5,...0.5]

46    continue

      f=0.d0
      do i=1,n-1
        f = f + 100.d0*(x(i+1)-x(i)+1.d0-x(i)*x(i))**2
      end do
*--
      g(1) = 200.d0*(x(2)-x(1)+1.d0-x(1)*x(1))*(-1.d0-2.d0*x(1))
      
      do i=2,n-1
        g(i) = 200.d0*(x(i)-x(i-1)+1.d0-x(i-1)*x(i-1))+
     *         200.d0*(x(i+1)-x(i)+1.d0-x(i)*x(i))*(-1.d0-2.d0*x(i))
      end do
      
      g(n) = 200.d0*(x(n)-x(n-1)+1.d0-x(n-1)**2)
      
      return





cF47                         ENGVAL1 (CUTE)                            
*                             
*                    Initial point x0=[2.,2.,2...,2.].
*
47    continue                          
      
      do i=1,n-1
        t(i) = x(i)*x(i) + x(i+1)*x(i+1)
      end do
      
      f = 0.d0
      
      do i=1,n-1
        f = f + t(i)*t(i) + (-4.d0*x(i) + 3.d0)
      end do   
*--
      g(1) = 4.d0*x(1)*t(1) - 4.d0
      
      do i=2,n-1
        g(i) = 4.d0*x(i)*t(i-1) + 4.d0*x(i)*t(i) - 4.d0
      end do
      
      g(n) = 4.d0*x(n)*t(n-1)  
      
      return




cF48                         DENSCHNA (CUTE)
*
*                       Initial point: [8, 8,...,8]
*
48    continue

      f=0.d0    
      do i=1,n/2
        f = f + x(2*i-1)**4 + 
     *          (x(2*i-1)+x(2*i))**2 + 
     *          (-1.d0+dexp(x(2*i)))**2
      end do
*--
      j=1      
      do i=1,n/2
        g(j) = 4.d0*x(2*i-1)**3 + 2.d0*(x(2*i-1)+x(2*i))
        g(j+1) = 2.d0*(x(2*i-1)+x(2*i)) + 
     *           2.d0*(dexp(x(2*i)))*(-1.d0+dexp(x(2*i)))    
        j=j+2
      end do
      
      return






cF49                   DENSCHNB  (CUTE)
*
*               Initial point: [0.1, 0.1,...,0.1]

49    continue

      f=0.d0
      do i=1,n/2
        f = f + (x(2*i-1)-2.d0)**2 +
     *          ((x(2*i-1)-2.d0)**2)*(x(2*i)**2) +
     *          (x(2*i)+1.d0)**2
      end do   
*--       
      j=1
      do i=1,n/2
        g(j) = 2.d0*(x(2*i-1)-2.d0) + 2.d0*(x(2*i-1)-2.d0)*x(2*i)*x(2*i)
        g(j+1) = ((x(2*i-1)-2.d0)**2)*2.d0*x(2*i) + 2.d0*(x(2*i)+1.d0)
        j=j+2
      end do
      
      return 




cF50                   DENSCHNC  (CUTE)
*
*               Initial point: [8, 8,...,8]

50    continue

      f=0.d0
      do i=1,n/2
        f = f + (-2.d0+x(2*i-1)**2+x(2*i)**2)**2 +
     *          (-2.d0+dexp(x(2*i-1)-1.d0)+x(2*i)**3)**2   
      end do
*--               
      j=1
      do i=1,n/2
       g(j) = 4.d0*x(2*i-1)*(-2.d0+x(2*i-1)**2+x(2*i)**2)+ 
     *  2.d0*(dexp(x(2*i-1)-1.d0))*(-2.d0+dexp(x(2*i-1)-1.d0)+x(2*i)**3)
       g(j+1) = 4.d0*x(2*i)*(-2.d0+x(2*i-1)**2+x(2*i)**2)+ 
     *  6.d0*(x(2*i)**2)*(-2.d0+dexp(x(2*i-1)-1.d0)+x(2*i)**3)
        j=j+2
      end do
      
      return





cF51                  DENSCHNF  (CUTE)
*
*               Initial point: [2,0,2,0,...,2,0]

51    continue

      f=0.d0
      do  i=1,n/2
        f=f+(2.d0*(x(2*i-1)+x(2*i))**2+(x(2*i-1)-x(2*i))**2-8.d0)**2+
     *      (5.d0*x(2*i-1)**2+(x(2*i)-3.d0)**2-9.d0)**2
      end do
*--
      j=1
      do i=1,n/2
      g(j)=2.d0*(2.d0*(x(2*i-1)+x(2*i))**2+(x(2*i-1)-x(2*i))**2-8.d0)*
     *  (4.d0*(x(2*i-1)+x(2*i))+2.d0*(x(2*i-1)-x(2*i))) +    
     *  2.d0*(5.d0*x(2*i-1)**2+(x(2*i)-3.d0)**2-9.d0)*10.d0*x(2*i-1)
      g(j+1)=2.d0*(2.d0*(x(2*i-1)+x(2*i))**2+(x(2*i-1)-x(2*i))**2-8.d0)*
     *  (4.d0*(x(2*i-1)+x(2*i))-2.d0*(x(2*i-1)-x(2*i))) +
     *  2.d0*(5.d0*x(2*i-1)**2+(x(2*i)-3.d0)**2-9.d0)*2.d0*(x(2*i)-3.d0)      
        j=j+2
      end do           
      return





cF52                     SINQUAD (CUTE)
*                    
*             Initial Point: [0.1, 0.1, ..., 0.1]

52    continue

      f=(x(1)-1.d0)**4 + (x(n)**2-x(1)**2)**2
      
      do i=1,n-2
        t(i) = dsin(x(i+1)-x(n)) - x(1)**2 + x(i+1)**2
        f = f + t(i)*t(i)
      end do
*--      
      g(1) = 4.d0*(x(1)-1.d0)**3 - 4.d0*x(1)*(x(n)**2-x(1)**2) 
      do i=1,n-2
        g(1) = g(1) - 4.d0*t(i)*x(1)
      end do
      
      do i=2,n-1
        g(i) = 2.d0*t(i-1)*(dcos(x(i)-x(n))+2.d0*x(i))
      end do
      
      g(n) = 4.d0*x(n)*(x(n)**2-x(1)**2)
      do i=1,n-2
        g(n) = g(n) - 2.d0*t(i)*dcos(x(i+1)-x(n))    
      end do

      return





cF53                         HIMMELBG (CUTE)
*
*                    Initial Point: [1.5,1.5,...,1.5]                 

53    continue

      f=0.d0
      do i=1,n/2
        f = f + (2.d0*x(2*i-1)**2+3.d0*x(2*i)**2)*
     *          (dexp(-x(2*i-1)-x(2*i)))
      end do
*--
      j=1
      do i=1,n/2
        t1 = 2.d0*x(2*i-1)**2+3.d0*x(2*i)**2
        t2 = dexp(-x(2*i-1)-x(2*i))          
        
        g(j)   =  4.d0*x(2*i-1)*t2 - t1*t2
                     
        g(j+1) =  6.d0*x(2*i)*t2   - t1*t2
                  
        j=j+2
      end do  
      return 




cF54                         HIMMELBH (CUTE)
*
*                    Initial Point: [1.5,1.5,...,1.5]                 

54    continue

      f= 0.d0
      do i=1,n/2
        f=f+(-3.d0*x(2*i-1)-2.d0*x(2*i)+2.d0+x(2*i-1)**3 + x(2*i)**2)
      end do
*--
      j=1
      do i=1,n/2
        
        g(j)   =  -3.d0 + 3.d0*x(2*i-1)**2
                     
        g(j+1) =  -2.d0 + 2.d0*x(2*i)
                  
        j=j+2
      end do  
      return  




cF55                       DIXON3DQ  (CUTE)
*
*                  Initial Point x0=[-1, -1,..., -1]
*
55    continue

      f=(x(1)-2.d0)**2
      
      do i=1,n-1
        f = f + (x(i)-x(i+1))**2
      end do
      
      f = f + (x(n)-1.d0)**2
*--
      g(1) = 2.d0*(x(1)-2.d0) + 2.d0*(x(1)-x(2))
      
      do i=2,n-1
        g(i) = -2.d0*(x(i-1)-x(i)) + 2.d0*(x(i)-x(i+1))
      end do
      
      g(n) = -2.d0*(x(n-1)-x(n)) + 2.d0*(x(n)-1.d0)
                
      return



cF56                     BIGGSB1 (CUTE)
*                    
*               Initial Point: [0., 0., ....,0.]

56    continue
      
      f=(x(1)-1.d0)**2 + (1.d0-x(n))**2
      do i=2,n
        f = f + (x(i)-x(i-1))**2
      end do
*--
      g(1) = 4.d0*x(1) - 2.d0*x(2) - 2.d0 
      
      do i=2,n-1
        g(i) = 4.d0*x(i) - 2.d0*x(i-1) - 2.d0*x(i+1)
      end do
      
      g(n) = 4.d0*x(n) - 2.d0*x(n-1) - 2.d0          
      
      return




cF57                       Perturbed Quadratic function
*
*                     Initial Point:  [0.5, 0.5, ......, 0.5].
*     
57    continue
      
      temp1 = 0.d0
      do i=1,n
        temp1 = temp1 + float(i)*x(i)
      end do  
      
      f = temp1*temp1
      
      do i=1,n
        f = f + float(i)*x(i)**2      
      end do  
*--      
      do i=1,n
        g(i) = float(i)*2.d0*x(i) + 2.d0*temp1*float(i)
      end do        
      return



cF58                          GENROSNB (CUTE)
*
*                    Initial Point: [-1.2, 1, ...  -1.2, 1]

58    continue

      f = (x(1)-1.d0)**2
      do i=2,n
        f = f + 100.d0*(x(i)-x(i-1)**2)**2
      end do
*--
      g(1) = 2.d0*(x(1)-1.d0)-400.d0*x(1)*(x(2)-x(1)**2)
      do i=2,n-1
        g(i) = 200.d0*(x(i)-x(i-1)**2)-400.d0*x(i)*(x(i+1)-x(i)**2)
      end do

      g(n) = 200.d0*(x(n)-x(n-1)**2)
            
      return




cF59                   QP1 Extended Quadratic Penalty
*
*                     Initial Point:  [1, 1, ......,1].
*                     

59    continue
      
      t1=0.d0
      do i=1,n
        t1 = t1 + x(i)*x(i)
      end do
      t1 = t1 - 0.5d0  
        
	f = 0.d0                                
      do i=1,n-1
        f = f + (x(i)*x(i) - 2.d0)**2
        g(i) = 4.d0*(x(i)*x(i)-2.d0)*x(i) + 4.d0*t1*x(i)
      end do  
      
      f = f + t1*t1
      g(n) = 4.d0*t1*x(n)
      
      return




cF60                 QP2 Quadratic Penalty QP2 Function           
*
*                     Initial Point:  [1, 1, ......,1].
*                     

60    continue
      
      t1=0.d0
      do i=1,n
        t1 = t1 + x(i)*x(i)
      end do
      t1 = t1 - 100.d0  
                   
      f = 0.d0                   
      do i=1,n-1
        f = f + (x(i)*x(i) - dsin(x(i)))**2
        g(i) = 2.d0*(x(i)*x(i)-dsin(x(i)))*(2.d0*x(i)-dcos(x(i))) 
     *       + 4.d0*t1*x(i)
      end do  
      
      f = f + t1*t1
      g(n) = 4.d0*t1*x(n)
      
      return





cF61                           STAIRCASE S1                               
*
*                      Initial point x0=[1,1,...,1].
*

61    continue

	f=0.d0
      do i=1,n-1
        f = f + (x(i)+x(i+1)-float(i))**2
      end do
*--
	  g(1) = 2.d0*(x(1)+x(2)-1.d0)
      do i=2,n-1
        g(i) = 2.d0*(x(i-1)+x(i)-float(i-1)) +  
     *         2.d0*(x(i)+x(i+1)-float(i))
      end do
      g(n) = 2.d0*(x(n-1)+x(n)-float(n-1))                     	  
      
      return





cF62                           STAIRCASE S2                              
*
*                      Initial point x0=[1,1,...,1].
*

62    continue
	f = 0.d0
	do i=2,n
	  f = f + (x(i-1)+x(i)-float(i))**2
	end do
*--
	g(1)=2.d0*(x(1)+x(2)-2.d0)
	do i=2,n-1
	  g(i) = 2.d0*(x(i-1)+x(i)-float(i)) + 
     *         2.d0*(x(i)+x(i+1)-float(i+1))
	end do
	g(n) = 2.d0*(x(n-1)+x(n)-float(n))

	return




cF63                           STAIRCASE S3                              
*
*                      Initial point x0=[2,2,...,2].
*

63    continue
	f = 0.d0
	do i=2,n
	  f = f + (x(i-1)+x(i)+float(i))**2
	end do
*--
	g(1) = 2.d0*(x(1)+x(2)+2.d0)

	do i=2,n-1
	  g(i) = 2.d0*(x(i-1)+x(i)+float(i))+
     *         2.d0*(x(i)+x(i+1)+float(i+1))
	end do

	g(n) = 2.d0*(x(n-1)+x(n)+float(n))

	return




cF64                          Extended Trigonometric ET2
*
*                          Initial Point: [0.2, 0.2, ....,0.2].

64    continue      
      s= float(n)
      do i=1,n
        s = s - dsin(x(i))
      end do  
	 
      do i=1,n
        temp(i) = s + float(i)*(1.d0-dsin(x(i))) - dsin(x(i))
      end do 

      f = 0.d0
      do i=1,n
        f = f + temp(i)**2
      end do   
*--
      s=0.d0
      do i=1,n
        s = s + temp(i)
      end do
      
      do i=1,n
        g(i) = -2.d0*s*dcos(x(i)) + 
     +          2.d0*temp(i)*(-float(i)*dcos(x(i))-dcos(x(i)))  
      end do      
      return



cF65                          QP3 Extended Quadratic Penalty 
*
*                          Initial Point: [1., 1., ....,1.].

65    continue

      t1=0.d0
      do i=1,n
        t1 = t1 + x(i)*x(i)
      end do
      t1 = t1 - 0.25d0  
        
	f = t1*t1  
	                              
      do i=1,n-1
        f = f - (x(i)*x(i) - 1.d0)**2
        g(i) = -4.d0*(x(i)*x(i)-1.d0)*x(i) + 4.d0*t1*x(i)
      end do  
      
      g(n) = 4.d0*t1*x(n)
      
      return	





cF66                          EG3                             
*                             
*                  Initial point x0=[1,1,1...,1].
*
*
66    continue

      f=0.5d0*dcos(x(n)*x(n))

      do i=1,n-1
        f = f + dcos(x(1)+x(i)*x(i)-1.d0)
      end do
*--
      g(1)=-(1.d0+2.d0*x(1))*dsin(x(1)+x(1)*x(1)-1.d0)
      do i=2,n-1
        g(1) = g(1) - dsin(x(1)+x(i)*x(i)-1.d0)
      end do               
      
      do i=2,n-1
        g(i) = -2.d0*x(i)*dsin(x(1)+x(i)*x(i)-1.d0)
      end do
      
      g(n) = -x(n)*dsin(x(n)*x(n))

      return




cF67                          GENROSEN-2                             
*                             
*                  Initial point x0=[-1.2, 1, -1.2, 1,...,1].
*
*
67    continue
	c = 100.d0

	f = (x(1)-1.d0)**2
	do i=2,n
	  f = f + c*(x(i-1)**2-x(i))**2
	end do
*--
	g(1) = 2.d0*(x(1)-1.d0) + 4.d0*c*(x(1)**2-x(2))*x(1)

	do i=2,n-1
	  g(i)=-2.d0*c*(x(i-1)**2-x(i))+4.d0*c*(x(i)**2-x(i+1))*x(i)
	end do

	g(n) = -2.d0*c*(x(n-1)**2-x(n))

	return





cF68                          PRODsin (m=n-1)                            
*                             
*                Initial point x0=[5. 5, 5, 5,...,5].
*
*
68    continue

	m = n-1
	t1=0.d0
	t2=0.d0

	do i=1,m	
	  t1 = t1 + x(i)*x(i)
	end do

	do i=1,n
	  t2 = t2 + dsin(x(i))
	end do

	f = t1*t2
*--
	do i=1,m
	  g(i) = 2.d0*x(i)*t2 + t1*dcos(x(i))
	end do

	do i=m+1,n
	  g(i) = t1*dcos(x(i))
	end do

	return






cF69                          PROD1 (m=n)                            
*                             
*              Initial point x0=[1. 1, 1, 1,...,1].
*
*
69    continue

	m = n
	t1=0.d0
	t2=0.d0

	do i=1,m	
	  t1 = t1 + x(i)
	end do

	do i=1,n	
	  t2 = t2 + x(i)
	end do

	f = t1*t2
*--
	do i=1,m
	  g(i) = t1+t2
	end do

	do i=m+1,n
	  g(i) = t1
	end do

	return





cF70                          PRODcos (m=n-1)                            
*                             
*                Initial point x0=[5. 5, 5, 5,...,5].
*
*
70    continue

	m = n-1
	t1=0.d0
	t2=0.d0

	do i=1,m	
	  t1 = t1 + x(i)*x(i)
	end do

	do i=1,n
	  t2 = t2 + dcos(x(i))
	end do

	f = t1*t2
*--
	do i=1,m
	  g(i) = 2.d0*x(i)*t2 - t1*dsin(x(i))
	end do

	do i=m+1,n
	  g(i) = -t1*dsin(x(i))
	end do

	return






cF71                          PROD2 (m=1)                            
*                             
*              Initial point x0=[15. 15, 15, 15,...,15].
*
*
71    continue

	m = 1
	t1=0.d0
	t2=0.d0

	do i=1,m	
	  t1 = t1 + x(i)**4
	end do

	do i=1,n	
	  t2 = t2 + float(i)*x(i)
	end do

	f = t1*t2
*--

	do i=1,m
	  g(i) = 4.d0*t2*x(i)**3 + float(i)*t1
	end do

	do i=m+1,n
	  g(i) = float(i)*t1
	end do

	return




cF72                          ARGLINB (m=5)                            
*                             
*         Initial point x0=[0.01   0.001, ....  ,0.01   0.001].
*
*
72    continue

      m=5
      f=0.d0         
      
      do i=1,m
        u(i)=0.d0
        do j=1,n
          u(i) = u(i) + float(i)*float(j)*x(j)
        end do  

        f = f + (u(i)-1.d0)**2
      end do  
*--
      do j=1,n
        g(j) = 0.d0
        do i=1,m
          g(j) = g(j) + 2.d0*(u(i)-1.d0)*float(i)*float(j)
        end do
      end do
                              
      return 





cF73                        DIXMAANA (CUTE)                             
*                             
*                    Initial point x0=[2.,2.,2...,2.].
*                           Modified m=n/4
73    continue
*    
        alpha = 1.d0
        beta  = 0.d0
        gamma = 0.125d0
        delta = 0.125d0
                      
        k1 = 0
        k2 = 0
        k3 = 0
        k4 = 0
           
        m = n/4

	f = 1.d0

	do i=1,n
	  f = f + alpha * x(i)*x(i)*((float(i)/float(n))**k1)
	end do

      do i=1,n-1
        f = f + beta * x(i)*x(i)*((x(i+1)+x(i+1)*x(i+1))**2) *
     *          ((float(i)/float(n))**k2)
      end do
        
      do i=1,2*m
        f = f + gamma * x(i)*x(i) * (x(i+m)**4) *
     *          ((float(i)/float(n))**k3)
      end do
        
      do i=1,m
        f = f + delta * x(i) * x(i+2*m) *
     *          ((float(i)/float(n))**k4)
      end do
*--
      do i=1,n
        g(i) =0.d0
      end do
c1        

      do i=1,n
        g(i) = g(i) + 2.d0*alpha*x(i)*((float(i)/float(n))**k1)  
      end do
c2	

	g(1) = g(1) + 2.d0*beta*x(1)*((x(2)+x(2)*x(2))**2)*
     *              ((float(1)/float(n))**k2)
      do i=2,n-1
        g(i) = g(i) + 2.d0*beta*(x(i-1)**2)*(x(i)+x(i)**2)*
     *                     (1.d0+2.d0*x(i))*((float(i-1)/float(n))**k2)+
     *                2.d0*beta*x(i)*((x(i+1)+x(i+1)**2)**2)*
     *                     ((float(i)/float(n))**k2)  
      end do
	g(n) = g(n) + 2.d0*beta*(x(n-1)**2)*(x(n)+x(n)**2)*
     *                     (1.d0+2.d0*x(n)) 
c3
      
      do i=1,2*m
        g(i) = g(i) + 2.d0*gamma*x(i)*(x(i+m)**4)* 
     *                ((float(i)/float(n))**k3)      
        g(i+m) = g(i+m) + gamma*(x(i)**2)*4.d0*(x(i+m)**3)*
     *                ((float(i)/float(n))**k3)   
      end do
c4
      
      do i=1,m
        g(i) = g(i) + delta*x(i+2*m)*((float(i)/float(n))**k4)
        g(i+2*m) = g(i+2*m) + delta*x(i)*((float(i)/float(n))**k4)
      end do  
          
      return		
                                      



cF74                        DIXMAANB (CUTE)                             
*                             
*                    Initial point x0=[2.,2.,2...,2.].
*                           Modified m=n/4
74    continue
*    
        alpha = 1.d0
        beta  = 0.0625d0
        gamma = 0.0625d0
        delta = 0.0625d0
                      
        k1 = 0
        k2 = 0
        k3 = 0
        k4 = 1
           
        m = n/4

	f = 1.d0

	do i=1,n
	  f = f + alpha * x(i)*x(i)*((float(i)/float(n))**k1)
	end do

      do i=1,n-1
        f = f + beta * x(i)*x(i)*((x(i+1)+x(i+1)*x(i+1))**2) *
     *          ((float(i)/float(n))**k2)
      end do
        
      do i=1,2*m
        f = f + gamma * x(i)*x(i) * (x(i+m)**4) *
     *          ((float(i)/float(n))**k3)
      end do
        
      do i=1,m
        f = f + delta * x(i) * x(i+2*m) *
     *          ((float(i)/float(n))**k4)
      end do
*--
      do i=1,n
        g(i) =0.d0
      end do
c1        

      do i=1,n
        g(i) = g(i) + 2.d0*alpha*x(i)*((float(i)/float(n))**k1)  
      end do
c2	

	g(1) = g(1) + 2.d0*beta*x(1)*((x(2)+x(2)*x(2))**2)*
     *              ((float(1)/float(n))**k2)
      do i=2,n-1
        g(i) = g(i) + 2.d0*beta*(x(i-1)**2)*(x(i)+x(i)**2)*
     *                     (1.d0+2.d0*x(i))*((float(i-1)/float(n))**k2)+
     *                2.d0*beta*x(i)*((x(i+1)+x(i+1)**2)**2)*
     *                     ((float(i)/float(n))**k2)  
      end do
	g(n) = g(n) + 2.d0*beta*(x(n-1)**2)*(x(n)+x(n)**2)*
     *                     (1.d0+2.d0*x(n)) 
c3
      
      do i=1,2*m
        g(i) = g(i) + 2.d0*gamma*x(i)*(x(i+m)**4)* 
     *                ((float(i)/float(n))**k3)      
        g(i+m) = g(i+m) + gamma*(x(i)**2)*4.d0*(x(i+m)**3)*
     *                ((float(i)/float(n))**k3)   
      end do
c4
      
      do i=1,m
        g(i) = g(i) + delta*x(i+2*m)*((float(i)/float(n))**k4)
        g(i+2*m) = g(i+2*m) + delta*x(i)*((float(i)/float(n))**k4)
      end do  
          
      return		




cF75                        DIXMAANC (CUTE)                             
*                             
*                    Initial point x0=[2.,2.,2...,2.].
*                           Modified m=n/4
75    continue
*    
        alpha = 1.d0
        beta  = 0.125d0
        gamma = 0.125d0
        delta = 0.125d0
                      
        k1 = 0
        k2 = 0
        k3 = 0
        k4 = 0
           
        m = n/4

	f = 1.d0

	do i=1,n
	  f = f + alpha * x(i)*x(i)*((float(i)/float(n))**k1)
	end do

      do i=1,n-1
        f = f + beta * x(i)*x(i)*((x(i+1)+x(i+1)*x(i+1))**2) *
     *          ((float(i)/float(n))**k2)
      end do
        
      do i=1,2*m
        f = f + gamma * x(i)*x(i) * (x(i+m)**4) *
     *          ((float(i)/float(n))**k3)
      end do
        
      do i=1,m
        f = f + delta * x(i) * x(i+2*m) *
     *          ((float(i)/float(n))**k4)
      end do
*--
      do i=1,n
        g(i) =0.d0
      end do
c1        

      do i=1,n
        g(i) = g(i) + 2.d0*alpha*x(i)*((float(i)/float(n))**k1)  
      end do
c2	

	g(1) = g(1) + 2.d0*beta*x(1)*((x(2)+x(2)*x(2))**2)*
     *              ((float(1)/float(n))**k2)
      do i=2,n-1
        g(i) = g(i) + 2.d0*beta*(x(i-1)**2)*(x(i)+x(i)**2)*
     *                     (1.d0+2.d0*x(i))*((float(i-1)/float(n))**k2)+
     *                2.d0*beta*x(i)*((x(i+1)+x(i+1)**2)**2)*
     *                     ((float(i)/float(n))**k2)  
      end do
	g(n) = g(n) + 2.d0*beta*(x(n-1)**2)*(x(n)+x(n)**2)*
     *                     (1.d0+2.d0*x(n)) 
c3
      
      do i=1,2*m
        g(i) = g(i) + 2.d0*gamma*x(i)*(x(i+m)**4)* 
     *                ((float(i)/float(n))**k3)      
        g(i+m) = g(i+m) + gamma*(x(i)**2)*4.d0*(x(i+m)**3)*
     *                ((float(i)/float(n))**k3)   
      end do
c4
      
      do i=1,m
        g(i) = g(i) + delta*x(i+2*m)*((float(i)/float(n))**k4)
        g(i+2*m) = g(i+2*m) + delta*x(i)*((float(i)/float(n))**k4)
      end do  
          
      return		





cF76                        DIXMAAND (CUTE)                             
*                             
*                    Initial point x0=[2.,2.,2...,2.].
*                           Modified m=n/4
76    continue
*    
        alpha = 1.d0
        beta  = 0.26d0
        gamma = 0.26d0
        delta = 0.26d0
                      
        k1 = 0
        k2 = 0
        k3 = 0
        k4 = 0
           
        m = n/4

	f = 1.d0

	do i=1,n
	  f = f + alpha * x(i)*x(i)*((float(i)/float(n))**k1)
	end do

      do i=1,n-1
        f = f + beta * x(i)*x(i)*((x(i+1)+x(i+1)*x(i+1))**2) *
     *          ((float(i)/float(n))**k2)
      end do
        
      do i=1,2*m
        f = f + gamma * x(i)*x(i) * (x(i+m)**4) *
     *          ((float(i)/float(n))**k3)
      end do
        
      do i=1,m
        f = f + delta * x(i) * x(i+2*m) *
     *          ((float(i)/float(n))**k4)
      end do
*--
      do i=1,n
        g(i) =0.d0
      end do
c1        

      do i=1,n
        g(i) = g(i) + 2.d0*alpha*x(i)*((float(i)/float(n))**k1)  
      end do
c2	

	g(1) = g(1) + 2.d0*beta*x(1)*((x(2)+x(2)*x(2))**2)*
     *              ((float(1)/float(n))**k2)
      do i=2,n-1
        g(i) = g(i) + 2.d0*beta*(x(i-1)**2)*(x(i)+x(i)**2)*
     *                     (1.d0+2.d0*x(i))*((float(i-1)/float(n))**k2)+
     *                2.d0*beta*x(i)*((x(i+1)+x(i+1)**2)**2)*
     *                     ((float(i)/float(n))**k2)  
      end do
	g(n) = g(n) + 2.d0*beta*(x(n-1)**2)*(x(n)+x(n)**2)*
     *                     (1.d0+2.d0*x(n)) 
c3
      
      do i=1,2*m
        g(i) = g(i) + 2.d0*gamma*x(i)*(x(i+m)**4)* 
     *                ((float(i)/float(n))**k3)      
        g(i+m) = g(i+m) + gamma*(x(i)**2)*4.d0*(x(i+m)**3)*
     *                ((float(i)/float(n))**k3)   
      end do
c4
      
      do i=1,m
        g(i) = g(i) + delta*x(i+2*m)*((float(i)/float(n))**k4)
        g(i+2*m) = g(i+2*m) + delta*x(i)*((float(i)/float(n))**k4)
      end do  
          
      return	





cF77                        DIXMAANL (CUTE)                             
*                             
*                    Initial point x0=[2.,2.,2...,2.].
*                           Modified m=n/4
77    continue
*    
        alpha = 1.d0
        beta  = 0.26d0
        gamma = 0.26d0
        delta = 0.26d0
                      
        k1 = 2
        k2 = 0
        k3 = 0
        k4 = 2
           
        m = n/4

	f = 1.d0

	do i=1,n
	  f = f + alpha * x(i)*x(i)*((float(i)/float(n))**k1)
	end do

      do i=1,n-1
        f = f + beta * x(i)*x(i)*((x(i+1)+x(i+1)*x(i+1))**2) *
     *          ((float(i)/float(n))**k2)
      end do
        
      do i=1,2*m
        f = f + gamma * x(i)*x(i) * (x(i+m)**4) *
     *          ((float(i)/float(n))**k3)
      end do
        
      do i=1,m
        f = f + delta * x(i) * x(i+2*m) *
     *          ((float(i)/float(n))**k4)
      end do
*--
      do i=1,n
        g(i) =0.d0
      end do
c1        

      do i=1,n
        g(i) = g(i) + 2.d0*alpha*x(i)*((float(i)/float(n))**k1)  
      end do
c2	

	g(1) = g(1) + 2.d0*beta*x(1)*((x(2)+x(2)*x(2))**2)*
     *              ((float(1)/float(n))**k2)
      do i=2,n-1
        g(i) = g(i) + 2.d0*beta*(x(i-1)**2)*(x(i)+x(i)**2)*
     *                     (1.d0+2.d0*x(i))*((float(i-1)/float(n))**k2)+
     *                2.d0*beta*x(i)*((x(i+1)+x(i+1)**2)**2)*
     *                     ((float(i)/float(n))**k2)  
      end do
	g(n) = g(n) + 2.d0*beta*(x(n-1)**2)*(x(n)+x(n)**2)*
     *                     (1.d0+2.d0*x(n)) 
c3
      
      do i=1,2*m
        g(i) = g(i) + 2.d0*gamma*x(i)*(x(i+m)**4)* 
     *                ((float(i)/float(n))**k3)      
        g(i+m) = g(i+m) + gamma*(x(i)**2)*4.d0*(x(i+m)**3)*
     *                ((float(i)/float(n))**k3)   
      end do
c4
      
      do i=1,m
        g(i) = g(i) + delta*x(i+2*m)*((float(i)/float(n))**k4)
        g(i+2*m) = g(i+2*m) + delta*x(i)*((float(i)/float(n))**k4)
      end do  
          
      return	





cF78                        VARDIM (CUTE)                             
*                             
*                    Initial point x0=[1-1/n, 1-2/n,...,1-n/n.].
*                           Modified m=n/4
78    continue

	s = float(n)*float(n+1)/2.d0

	t1=0.d0
	do i=1,n
	  t1 = t1 + float(i)*x(i)
	end do
	t1 = t1-s

	f = 0.d0
	do i=1,n
	  f = f + (x(i)-1.d0)**2
	end do

	f = f + t1**2 + t1**4
*--
	do i=1,n
	  g(i) = 2.d0*(x(i)-1.d0) + 2.d0*t1*float(i) + 4.d0*float(i)*t1**3
	end do

	return





cF79                             DIAG-AUP1                             
*                             
*                    Initial point x0=[4., 4., ....4.].
*

79    continue                  

      f=0.d0
      do i=1,n
        f = f + 4.d0*(x(i)*x(i) - x(1))**2 + (x(i)**2-1.d0)**2
      end do  

*--
      g(1) = 4.d0*(x(1)**2-1.d0)*x(1) + 
     *       8.d0*(x(1)*x(1)-x(1))*(2.d0*x(1)-1.d0)
      do i=2,n
        g(1) = g(1) - 8.d0*(x(i)*x(i)-x(1))
      end do
      
      do i=2,n
        g(i) = 16.d0*x(i)*(x(i)*x(i)-x(1)) + 4.d0*(x(i)**2-1.d0)*x(i)
      end do  
        
      return





cF80                             ENGVAL8                             
*                             
*                    Initial point x0=[2., 2., ....2.].
*
80	continue
	
	f = 0.d0
	do i=1,n-1
	  f = f + (x(i)**2+x(i+1)**2)**2 - (7.d0-8.d0*x(i))
	end do
*--
	g(1) = 4.d0*(x(1)**2+x(2)**2)*x(1) + 8.d0

	do i=2,n-1
	  g(i)=4.d0*(x(i-1)**2+x(i)**2)*x(i) + 
     *       4.d0*(x(i)**2+x(i+1)**2)*x(i) + 8.d0
	end do

	g(n) = 4.d0*(x(n-1)**2+x(n)**2)*x(n)

	return





***                  
      end

c------------------------------------------------ End EVALFG
*  
