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

      program graspmain

C     PARAMETERS
      integer nmax,mmax, maxiter
      parameter ( nmax    = 1000 )
      parameter ( mmax    =    1 )
      parameter ( maxiter = 2000 )

C     COMMON SCALARS

C     COMMON ARRAYS

C     LOCAL SCALARS
      real ntime
      integer inform,iter,n,m,i,itbest
      double precision gpsupn,seed,drand,tmax,f,vdummy
C    + drand

C     LOCAL ARRAYS
      real dum(2)
      double precision g(nmax),l(nmax),u(nmax),x(nmax),
     + lambda(mmax),rho(mmax),equatn(mmax),linear(mmax)

C     COMMON BLOCKS

C     SET PROBLEM DATA

C     SET LOCAL SCALARS
C      call initproblem(n,m,dlx,dux)
      call inip(n,x,l,u,m,lambda,rho,equatn,linear)

C     SET MAXIMUM CPU TIME

C      write(*,*) 'Enter the maximum CPU time (in seconds): '
C      read(*,*) tmax
      tmax = 100000

C     MAIN LOOP: TRYING TO PACK AS MANY AS POSSIBLE

      iter        =   0

C     SET A RANDOM INITIAL GUESS

C     Set the seed for the initial guess random generation

      seed = iter

C     Random generation of the initial guess

      do i = 1,n
C          x(i) = l(i) + ( u(i) - l(i) ) * drand(seed)
          x(i) = l(i) + ( u(i) - l(i) ) * 0.5d0
C          x(i) = 6.5d0
      end do

C      call evalal(n,x,m,f,inform)
      call evalf(n,x,f,inform)
      call evalgsupn(n,x,g,l,u,gpsupn)

      write(*,*) '--------------- INICIAL -------------------'
      write(*,*) "t | fg > ", ntime,' | ',
     +            f,gpsupn
      write(*,*) "| x | "
      do i = 1,n
          write(*,*) x(i)
      end do
 
 100  continue

C     TEST STOOPING CRITERIA

      if ( ntime .gt. tmax ) then
          write(*,*) 'GRASP Flag: CPU Time is exhausted.'
          go to 500
      end if

      iter = iter + 1

C     CALL THE OPTIMIZER

      call tun_grasp(n,x,m,l,u,f,g,gpsupn,iter,inform,itbest,ntime)

C     ITERATE

      if ( iter .gt. maxiter ) then
          write(*,*) 'GRASP Flag: Max trial is exhausted.'
          go to 500
      else
          seed = iter+1
          vdummy = drand(seed)

          do i = 1,n
              x(i) = l(i) + drand(seed) * ( l(i) + u(i) ) / 2.0d0
          end do

      end if

      go to 100

C     END OF MAIN LOOP

 500  continue

      call evalgsupn(n,x,g,l,u,gpsupn)

      ntime = etime(dum)

      write(*,*) '--------------- FINAL -------------------'
      write(*,*) "t | it | fg > ", ntime,' | ',itbest,' | ',
     +            f,gpsupn
      write(*,0010) f
      write(*,*) "| x | "
      do i = 1,n
          write(*,*) x(i)
      end do

      stop

C     NON-EXECUTABLE STATEMENTS
 0010 format('f = ',D25.16)
 0011 format(I7,1X,D25.16,1X,D25.16,1X,F8.2)

      end

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

      double precision function drand(ix)
      double precision ix

      double precision a,p,b15,b16,xhi,xalo,leftlo,fhi,k
      data a/16807.d0/,b15/32768.d0/,b16/65536.d0/,p/2147483647.d0/

      xhi= ix/b16
      xhi= xhi - dmod(xhi,1.d0)
      xalo= (ix-xhi*b16)*a
      leftlo= xalo/b16
      leftlo= leftlo - dmod(leftlo,1.d0)
      fhi= xhi*a + leftlo
      k= fhi/b15
      k= k - dmod(k,1.d0)
      ix= (((xalo-leftlo*b16)-p)+(fhi-k*b15)*b16)+k
      if (ix.lt.0) ix= ix + p
      drand= ix*4.656612875d-10

      return

      end

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

      subroutine evalgsupn(n,x,g,l,u,pgpsupn)

C     This routine was a copy from Genpack's:
C     "Compute continuous-project-gradient Euclidian and Sup norms,
C     internal gradient Euclidian norm, and store in nind the number of
C     free variables and in array ind their identifiers."

C     On Entry:

C     n      integer
C            number of variables
C
C     x      double precision x(n),
C            estimation of the solution,
C
C     g      double precision g(n)
C            gradient of the objective function at x
C
C     l      double precision l(n)
C            lower bounds on the variables
C
C     u      double precision u(n)
C            upper bounds on the variables
C
C
C     On Return:
C
C     pgpsupn   double precision
C              Squared Euclidian sup-norm of the continuous projected 
C              gradient at the final estimation

C     SCALAR ARGUMENTS
      integer n
      double precision pgpsupn 

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

C     LOCAL SCALARS
      integer k
      double precision pgpi

      pgpsupn = 0.0d0

      do k = 1,n
          pgpi     = min( u(k), max( l(k), x(k) - g(k) ) ) - x(k)
          pgpsupn  = max( pgpsupn, abs( pgpi ) )
      end do


      end