	PROGRAM GSA
!  * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * |
!     This routine start the GSA-loop

      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (MaxDim=500) 

      COMMON /Par1/ qV1,qV2,qA1,qT1,D,exp1,exp2,Coef,Tqt,To,ToScale
      COMMON /Par2/ qA,qV,qT
      COMMON /Par3/ NDimension,NStopMax,NRAN
      COMMON /Xvector/ X_t(Maxdim), X_0(Maxdim), X_Min(Maxdim)
      DATA One /1.0D+00/

	CALL GSAini()  ! Read initial GSA parameters in the file "GSA.in"

! IMPORTANT!  Is it necessary to initialize the vectors X_0 e X_Min.

	DO i=1,NDimension
	   X_Min(i) = X_0(i)
	   X_t(i)   = X_0(i)
	ENDDO
          func_0   = func(X_0)  
	func_Min = func_0
	func_t   = func_0

      OneqA1 = One/qA1
      Time    = 0.0D0
      NCycle  = 0

! Starting minimization loop

      DO WHILE (NCycle.LE.NStopMax)

         Time = Time + One 
         NCycle = NCycle + 1

         T = Tqt/((One+Time)**qT1 - One) ! Temperature T(t)
         IF(D.EQ.0.0D0) THEN
            Tup = One
         ELSE
            Tup = T**(D/(qT-3.0D0))
         ENDIF

!--Create the new vector  X_t  using the g(qV,t) function.
         CALL Gfunc(T,Tup)

!--Evaluate the Potential function in terms of X_t.
	   func_t = func(X_t)

!--Verify if the Potential(X_n+1) is going down. In this case
!  change  X(n) ---> X(n+1). In the other case retains X(n).

         IF(func_t .LE. func_0) THEN 
             DO I=1, NDimension  
                X_0(I) = X_t(I) 
             ENDDO  
             func_0  = func_t 
             IF(func_t .LE. func_Min) THEN 
                func_Min = func_t               
                DO I=1,NDimension
                   X_Min(I) = X_t(I) 
                ENDDO 
             ENDIF 
         ELSE 
!            Evaluate the Acceptance Probability Function [0,1]. 
             DeltaE = func_t - func_0
!---         Xiang (Phys. Let. A 233, (1997) 216-220)
!		   qA = -3 -0.85*Time
		   qA1= qA - One
!---
!            PqA = One/(One+(One+qA1*DeltaE/T)**OneqA1)
             PqA = One/((One+qA1*DeltaE/T)**OneqA1)  ! Correo
             Rand = RAN3(NRAN)
!            If  Rand > Acc_Prob  retains   X_0. 
!            If  Rand < Acc_Prob  replace   X_0 by the new X_t. 
             IF(Rand .LT. PqA) THEN 
                DO I=1,NDimension 
                   X_0(I) = X_t(I) 
                ENDDO 
                func_0 = func_t 
             ENDIF 
         ENDIF 
      ENDDO

      END

!  * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * | 
      SUBROUTINE Gfunc(T,Tup)
!  * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * |
!     This routine evaluate the new set o parameters using the  
!     visiting probability function g(q,T).

      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (MaxDim=500) 
      COMMON /Par1/ qV1,qV2,qA1,qT1,D,exp1,exp2,Coef,Tqt,To,ToScale
      COMMON /Par3/ NDimension,NStopMax,NRAN
      COMMON /Xvector/ X_t(Maxdim), X_0(Maxdim), X_Min(Maxdim)

      DATA One /1.00000D+00/

      DO I = 1,NDimension
 		R = RAN3(NRAN)
		S = RAN3(NRAN)
		DeltaX = Coef*Tup/(One+qV1*R*R/T**exp1)**exp2  
		IF(S.LE.0.5) DeltaX = -DeltaX
		X_t(I) = X_0(I) + DeltaX				
      ENDDO


      RETURN
      END

!  * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * |
      FUNCTION dgamma(r)
!  * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * |
!     This routine evaluate usual Gamma function. 

      IMPLICIT REAL*8 (a - h, o - z)

      PARAMETER (
     1  p0 =  0.999999999999999990d+00, p1 = -0.422784335098466784d+00, 
     2  p2 = -0.233093736421782878d+00, p3 =  0.191091101387638410d+00, 
     3  p4 = -0.024552490005641278d+00, p5 = -0.017645244547851414d+00, 
     4  p6 =  0.008023273027855346d+00)
      PARAMETER (
     1  p7  = -0.000804329819255744d+00,p8 =  -0.000360837876648255d+00, 
     2  p9  =  0.000145596568617526d+00,p10 = -0.000017545539395205d+00, 
     3  p11 = -0.000002591225267689d+00,p12 =  0.000001337767384067d+00, 
     4  p13 = -0.000000199542863674d+00)

      n = NINT(r - 2)
      w = r - (n + 2)
      y = ((((((((((((p13   * w + p12)* w + p11)* w + p10)* 
     1    w + p9) * w + p8) * w + p7) * w + p6) * w + p5) * 
     2    w + p4) * w + p3) * w + p2) * w + p1) * w + p0
      IF (n .gt. 0) THEN
          w = r - 1
          DO k = 2, n
              w = w * (r - k)
          END DO
      ELSE
          w = 1
          DO k = 0, -n - 1
              y = y * (r + k)
          END DO
      END IF
      dgamma = w / y

	RETURN
      END

!    * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * |
      DOUBLE PRECISION FUNCTION RAN3(IDUM)
!    * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * |
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
!      This routine return an randomic number   0<= r <= 1

!     IMPLICIT REAL*4(M)
!     PARAMETER (MBIG=4000000.,MSEED=1618033.,MZ=0.,FAC=2.5D-7)
      PARAMETER (MBIG=1000000000,MSEED=161803398,MZ=0,FAC=1.D-9)
      COMMON /RAN1A/ IFF,MJ,MK,INEXT,INEXTP,MA(55)
!      DATA IFF /0/
!           INEXT=0
!           INEXTP=31
 
      IF(IDUM.LT.0.OR.IFF.EQ.0)THEN
      IFF=1
      MJ=MSEED-IABS(IDUM)
      MJ=MOD(MJ,MBIG)
      MA(55)=MJ
      MK=1
      DO 11 I=1,54
         II=MOD(21*I,55)
         MA(II)=MK
         MK=MJ-MK
         IF(MK.LT.MZ)MK=MK+MBIG
         MJ=MA(II)
11    CONTINUE
      DO 13 K=1,4
         DO 12 I=1,55
            MA(I)=MA(I)-MA(1+MOD(I+30,55))
            IF(MA(I).LT.MZ)MA(I)=MA(I)+MBIG
12       CONTINUE
13    CONTINUE
      INEXT=0
      INEXTP=31
      IDUM=1
      ENDIF
 
      INEXT=INEXT+1
      IF(INEXT.EQ.56)INEXT=1
      INEXTP=INEXTP+1
      IF(INEXTP.EQ.56)INEXTP=1
      MJ=MA(INEXT)-MA(INEXTP)
      IF(MJ.LT.MZ)MJ=MJ+MBIG
      MA(INEXT)=MJ
      RAN3=MJ*FAC
 
      RETURN
      END

!  * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * |
      SUBROUTINE GSAini()
!  * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * |
!     This routine read and initialize the GSA parameters
!     The number of parameters to be optimized (NDIMENSION) is problem
!     dependent.

      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (MaxDim=500) 
      COMMON /Par1/ qV1,qV2,qA1,qT1,D,exp1,exp2,Coef,Tqt,To,ToScale
	COMMON /Par2/ qA,qV,qT
      COMMON /Par3/ NDimension,NStopMax,NRAN

      DATA One /1.00000D+00/
	CHARACTER tempo*8, lixo*10 

	OPEN(7,FILE="GSA.in")
          READ(7,"(A10)") lixo ! Skip the comment first line
		READ(7,*) qA         ! Acceptance index 
		READ(7,*) qT		 ! Temperature index
		READ(7,*) qV		 ! Visiting index
		READ(7,*) NStopMax   ! Max number of GSA-loops	
		READ(7,*) To         ! Initial Temperature
		READ(7,*) nDimension ! Number of parameters to be optimized
	CLOSE(7)

	CALL TIME(tempo)

	READ(tempo(7:8),"(I2)") NRAN   ! Seed time dependent
	NRAN = -NRAN

	D = 0  ! Dimension problem 

      Pi = 3.14159265359D0

!    Acceptance probability
      qa1= qa - one

!    Temperature
      qT1 = qT - One
      Tqt = To*(2.0D0**qT1-One)

!    Visiting probability
      qv1= qv - One
      qv2= 2.0D0**qv1 - One
      tmp= One/qv1 - 0.5D0
      GamaDown = dgamma(tmp)
      exp1 = 2.0D0/(3.0D0 - qv)

      IF(D.EQ.0.0D0) THEN
         Coef1  = One 
         exp2   = One/qv1 - 0.5D0
         GamaUp = GamaDown
      ELSE
         Coef1  = (qv1/Pi)**(D*0.5D0)
         exp2   = One/qv1 + 0.5D0*D - 0.5D0
         GamaUp = dgamma(exp2)
      ENDIF

      Coef = Coef1*GamaUp/GamaDown


      RETURN
      END


!    * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * |
      FUNCTION func(X)
!    * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * |	
	IMPLICIT DOUBLE PRECISION (A-H,O-Z)
	PARAMETER (MaxDim=500)
      COMMON /Par3/ NDimension,NStopMax,NRAN

      DIMENSION X(MaxDim)

!	CALL here your routine with the objective function


!     X(i) is parameters set and f is the value of the objective function
	CALL funct(X,f)
	func = f

      RETURN	
      END

