C     =================================================================
C     File: c-grasp.f
C     =================================================================

      program c_gencan

C     This program try to find some y such that f(y) <= f(x) by
C     using the C-GRASP algorithm, described in:
C
C     M.J. HIRSCH, P.M. PARDALOS, M.G.C. RESENDE. SPEEDING UP
C     CONTINUOUS GRASP. AT&T Labs Research Technical Report TD-6U2P2H, 
C     AT&T Shannon Laboratory, 180 Park Avenue, Florham Park, 
C     NJ 07932 USA, September 2006.
C
C     and
C
C     M.J. HIRSCH, C.N. Meneses, P.M. PARDALOS, M.G.C. RESENDE. Global 
C     optimization by continuous GRASP. AT&T Labs Research Technical 
C     Report TD-6MPUV9. AT&T Labs Research, Florham Park, NJ 07932, 
C     March 8, 2006.

C     PARAMETERS
      integer MaxIter,nmax,mmax
      double precision INFINITO,alpha
      parameter ( MaxIter             = 2000       )
      parameter ( INFINITO = 10000000000000000.0d0 )
      parameter ( alpha               = 0.3d0      )
      parameter ( nmax    = 500000 )
      parameter ( mmax    = 500000 )

C     LOCAL SCALARS
      real etime,ntime,btime,ctime
      integer n,m,i,iter,itbest,bfound,
     + nsamples,bnsamples,totevalf,bevalf,totevalg,bevalg,fcnmb
      double precision seed,drand,h,BestSolValue,SolValue,difSol,
     + bdifSol,vdummy,h_s,h_e
      logical impr_c,impr_l

C     LOCAL ARRAYS
      real dum(2)
      double precision BestSolution(nmax),y(nmax),
     + l(nmax),u(nmax)
 
C     GENCAN COMPATIBILITY

C     GENCAN'S SCALARS
      logical checkder
      character * 6 precond
      integer gtype,hptype,iprint,maxoutit,maxtotfc,maxtotit,ncomp,
     + inform,outiter,totcgcnt,totfcnt,totgcnt,totiter
      double precision epsfeas,epsopt,nalpsupn,snorm
      real time

C     GENCAN'S ARRAYS
      integer wi1(nmax)
      logical equatn(mmax),linear(mmax)
      double precision lambda(mmax),rho(mmax),wd1(mmax),
     +        wd2(mmax),wd3(nmax),wd4(mmax),wd5(mmax),wd6(nmax),
     +        wd7(nmax),wd8(nmax),wd9(nmax),wd10(nmax),wd11(nmax),
     +        wd12(nmax),wd13(nmax),wd14(nmax),wd15(nmax),wd16(nmax),
     +        wd17(nmax),wd18(nmax),wd19(nmax)

C     EXTERNAL SUBROUTINES
      external solver

C     SETUP PROBLEM

      call inip(n,y,l,u,m,lambda,rho,equatn,linear)

C     SET LOCAL SCALARS

      h_s = max(1.0d0, (u(1)-l(1)) / 1000.0d0 )
      h_e = 0.1d0 * h_s

C      write(*,*) '>>>> Inicio <<<< | hs,he = ',h_s,h_e

      write(*,*) " "
      write(*,*) " "
      write(*,0020) 
      write(*,0021) 
      write(*,0022) 

      SolValue     = INFINITO
      BestSolValue = INFINITO

C     SET LOCAL COUNTERS

      totevalf = 0
      totevalg = 0
      nsamples = 0
      bfound   = 0

C    SET TIME CHECK COUNTER

      ctime = etime(dum)

C     SET SOME SOLVER ARGUMENTS

      gtype   = 0
      hptype  = 6
      precond = 'NONE'
C      precond = 'QNCGNA'

      checkder = .false.

      epsfeas = 1.0d-04
      epsopt  = 1.0d-04

      maxoutit = 50
      maxtotit = 1000000
      maxtotfc = 5 * maxtotit

      iprint = 0
      ncomp  = 5

C     START METHOD

      iter         = 1

C     TEST STOP CRITERIA

      do while ( iter .le. MaxIter )

C         SET INITIAL RANDOM POINT

          seed = iter
          vdummy = drand(seed)

          do i = 1,n
               y(i) = l(i) + ( u(i) - l(i) ) * drand(seed)
          end do
          nsamples = nsamples + 1

          h = h_s

C         TEST GRID DENSITY

          do while( ( h .ge. h_e ) .and.
     +              ( iter .le. MaxIter ) )

              impr_c = .false.
              impr_l = .false.

C             C-GRASP-GENCAN GLOBAL PHASE

              call ConstructGreedyRandomizedSolution(y,n,h,l,u,seed,
     +         alpha,fcnmb,impr_c)

              call evalf(n,y,difSol,inform)

C             C-GRASP-GENCAN LOCAL PHASE - CALL OPTIMIZATION SOLVER

              call solver(n,y,l,u,m,lambda,rho,equatn,linear,gtype,
     + hptype,precond,checkder,epsfeas,epsopt,maxoutit,maxtotit,
     + maxtotfc,iprint,ncomp,SolValue,snorm,nalpsupn,outiter,totiter,
     + totfcnt,totgcnt,totcgcnt,
     + time,inform,wi1,wd1,wd2,wd3,wd4,wd5,wd6,wd7,wd8,wd9,wd10,wd11,
     + wd12,wd13,wd14,wd15,wd16,wd17,wd18,wd19)

              totevalf = totevalf + fcnmb + totfcnt
	      totevalg = totevalg + totgcnt

C             TEST SOLUTION IMPROVEMENT IN LOCAL PHASE

              difSol = SolValue - difSol

              if( abs(difSol) .gt. 0.0d0 ) then
                  impr_l = .true.
              end if

C             UPDATES THE BEST SOLUTION

              ntime = etime(dum)

	      if( BestSolValue .gt. 
     +            ( SolValue + 0.0001d0 * abs( SolValue )
     +            + 1.0E-6 ) ) then	  

                  do i = 1,n
                      BestSolution(i) = y(i)
                  end do

                  itbest    = iter
                  btime     = ntime
                  bevalf    = totevalf
	          bevalg    = totevalg
	          bnsamples = nsamples
	          bdifSol   = difSol
	          bfound    = bfound + 1

           write(*,0023) bfound,iter,SolValue,btime,bdifSol

                  BestSolValue = SolValue

              end if

C             CHECK TIME SPENT

              if( ntime - ctime .gt. 60.0d0 ) then

                  ctime = ntime
                  write(*,0010) ctime,iter,h

              end if

C             MAKES THE GRID MORE DENSE

              if( .not.( impr_c ) .and. .not.( impr_l ) ) then
                  h = h / 2.0d0
              end if

C             ITERATE

              iter = iter + 1

          end do

      end do

C     WRITES THE BEST SOLUTION FOUND

      write(*,*) " "
      write(*,*) " "
      write(*,*) " "
      write(*,0030)
      write(*,0031) BestSolValue,BestSolValue
      write(*,0035)
      write(*,0032) btime,ntime
      write(*,0033) itbest,bnsamples
      write(*,0034) bevalf,bevalg
C      write(*,0035)
C      write(*,*) " "
C      do i = 1,n
C         write(*,0036) i,BestSolution(i)
C      end do
      write(*,0037)

      write(*,*) " "
      write(*,*) '* Original experiments executed with G77. This ',
     + 'result may reach a diferent soluction since G77 and ',
     + 'gfortran are not fully compatible. *'

C     NON-EXECUTABLE STATEMENTS
 0010 format("---- TIME : ",F7.2," -- Iter : ",I4,
     + " -- h : ",F9.4," ----")
 0020 format("-----------------< Best solution found >",
     + "----------------")
 0021 format(" #    Iter                  f*     t*",
     + "      GENCAN Improv")
 0022 format("---  -----  ------------------  -----",
     + "  -----------------")
 0023 format(I3,2X,I5,2X,F18.10,2X,F4.2,2X,F18.10)
 0030 format("******************** FINAL SOLUTION ",
     + "********************")
 0031 format("f(x*) : ",D25.16," ( ",F18.10," )" )
 0032 format("time(x*) : ",F12.5," | total time : ",F12.5 )
 0033 format("iter(x*) : ",I4,8X," | #x0        : ",I7 )
 0034 format("eval_f   : ",I7,5X," | eval_g     : ",I7 )
 0035 format("------------------------------------",
     + "--------------------")
 0036 format("x( ",I4," ) = ",D25.16)
 0037 format("************************************",
     + "********************")

      end;

C     ******************************************************************
C     ******************************************************************

      subroutine ConstructGreedyRandomizedSolution(x,n,h,l,u,seed,
     +                                   alpha,totevalf,impr_c)
C     This subroutine try to find some y such that f(y) <= f(x) by
C     using the C-GRASP algorithm, described in:
C
C     M.J. Hirsch, P.M. Pardalos, and M.G.C. Resende. 
C     Speeding up continuous GRASP.
C     Technical Report TD-6U2P2H, AT&T Labs Research, Florham Park, NJ 07932,
C     2006.
C
C     M.J. Hirsch, C.N. Meneses, P.M. Pardalos, and M.G.C. Resende. 
C     Global optimization by continuous GRASP. 
C     Optimization Letters, 1:201–212, 2007.
C
C     C.N. MENESES, P.M. PARDALOS, AND M.G.C. RESENDE. GRASP FOR
C     NONLINEAR OPTIMIZATION. AT&T Labs Research Technical
C     Report TD-6DUTRG.
C
C     Parameters of the subroutine:
C
C     On Entry:
C
C     x        double precision x(n) : Initial point x^0
C
C     n        integer : number of variables
C
C     h        double precision : grid step length
C
C     l        integer : intervals lower bound
C
C     u        integer : intervals upper bound
C
C     seed     double precision : randomic factor (used by drand)
C
C     alpha    double precision : heuristic RCL constant
C
C     On Return:
C
C     x        double precision x(n) : min(x^0,y)
C
C     seed     double precision : randomic factor (used by drand)
C
C     totevalf integer : total functions evaluations
C
C     impr_c   boolean : true if y < x^0; false, otherwise
C
C     PARAMETERS
      double precision INFINITO
      parameter ( INFINITO = 10000000000000000.0d0 )

C     SCALAR ARGUMENTS
      integer n,totevalf,fcnmb
      double precision seed,alpha,h
      logical impr_c,ReUse

C     ARRAY ARGUMENTS
      double precision x(n),l(n),u(n)

C     LOCAL SCALARS
      integer i,j,k,RandomlySelectElement
      double precision vmax,vmin,xibest

C     LOCAL ARRAYS
      integer S(n+1),RCL(n+1)
      double precision z(n),g(n),y(n)

C     SET LOCAL SCALARS

      totevalf = 0
      fcnmb = 0

      do i = 1,n

          S(i)   = i
          RCL(i) = -1
          y(i)   = x(i)

      end do

      S(n+1)   = n
      RCL(n+1) = 0

C     TRY TO FIND SOME POINT z SUCH THAT f(x) < f(z)

      ReUse = .false.

      do while( S(n+1) .gt. 0 )

          vmin = INFINITO
          vmax = ( -1.0d0 ) * INFINITO

          do i = 1,S(n+1)

C             MINIMIZE COMPONENT y( S(i) ) WITH STEP h

              if( .not.( ReUse ) ) then

C                      Minimize(x,h,indx,n,l,u,xibest,fbest,m)
                  call Minimize(y,h,S(i),n,l,u,xibest,g(S(i)),fcnmb)
	          totevalf = totevalf + fcnmb

                  z(S(i)) = xibest

              end if

              if( vmin .gt. g(S(i)) ) then

                  vmin = g(S(i))

              end if

              if( vmax .lt. g(S(i)) ) then

                  vmax = g(S(i))

              end if

          end do

C         SET SOLUTION CANDIDATES

          RCL(n+1) = 0

          do i=1,S(n+1)

              if( g( S(i) ) .le. ( vmin + alpha * ( vmax - vmin ) ) )
     +        then

                  RCL(n+1) = RCL(n+1) + 1
                  RCL( RCL(n+1) ) = i

              end if

          end do

C         RANDOMLY SELECT A RCL ELEMENT

          seed = ( S(n+1) + seed ) / 2.0d0
          j = RandomlySelectElement(seed,1,RCL(n+1))

C         UPDATES x,y

          if( y( RCL(j) ) .eq. z( RCL(j) ) ) then
              ReUse = .true.
          else
              x( RCL(j) ) = z( RCL(j) )
              y( RCL(j) ) = z( RCL(j) )
              ReUse = .false.
              impr_c = .true.
          end if

C         S = S \ {j}

          k           = S( RCL(j) )
          S( RCL(j) ) = S( S(n+1) )
          S( S(n+1) ) = k
          S(n+1)      = S(n+1) - 1

      end do

      end

C     ******************************************************************
C     ******************************************************************

      integer function RandomlySelectElement(seed,l,u)
C     This subroutine returns a random integer between l and u

C     Parameters of the subroutine:
C
C     On Entry:
C
C     seed     double precision
C
C     l        integer : intervals lower bound
C
C     u        integer : intervals upper bound
C
C     On Return:
C
C     function integer : returns i \in [l,u]
C
C     seed     double precision
C
C     SCALAR ARGUMENTS
      integer l,u
      double precision seed

C     LOCAL SCALARS
      integer delta
      double precision drand,nran

      delta = u - l + 1
      nran  = drand(seed)
      nran  = drand(seed)

      RandomlySelectElement = max(nint( nran * delta ),1)

      return

      end;

C     ******************************************************************
C     ******************************************************************

      subroutine Minimize(x,h,indx,n,l,u,xibest,fbest,totevalf)
C
C     Parameters of the subroutine:
C
C     On Entry:
C
C     x        double precision x(n),
C              initial point,
C
C     h        double precision,
C              grid density
C
C     indx     integer,
C              index to be moved,
C
C     n        integer,
C              number of variables,
C
C     l        double precision l(n),
C              lower bounds on x,
C
C     u        double precision u(n),
C              upper bounds on x,
C
C     On Return
C
C     xibest   double precision,
C              best position found
C
C     fbest    double precision
C              best f's value (with xibest)

C     PARAMETERS
      double precision INFINITO
      parameter ( INFINITO = 10000000000000000.0d0 )

C     SCALAR ARGUMENTS
      integer n,totevalf,indx
      double precision h,xibest,fbest

C     ARRAY ARGUMENTS
      double precision x(n),l(n),u(n)

C     LOCAL SCALARS
      integer i,err
      double precision vmin,vle

C     LOCAL ARRAYS
      double precision t(n)

C     SET LOCAL SCALARS

      totevalf = 0

      do i = 1,n

          t(i) = x(i)

      end do

C     SET STARTUP MIN VALUE

      call evalf(n,t,vle,err)
      totevalf = totevalf + 1
      vmin = vle
      xibest = t(indx)

      t(indx) = l(indx)

C     SEARCH FOR A BEST VALUE

      do while ( t(indx) .le. u(indx) )

          call evalf(n,t,vle,err)
	  totevalf = totevalf + 1

          if( vmin .gt. vle ) then

              vmin = vle
              xibest = t(indx)

          end if

          t(indx) = t(indx) + h

      end do

      fbest = vmin

      end

