      SUBROUTINE PLUG1
C The Price algorithm as modified by F. Theos, D.G. Papageorgiou and I. E. Lagaris
      call prcglo
      END

      SUBROUTINE PRCGLO    
* This is an implementation of the improved Price algorithm                    
* for global optimization.
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)      
* MAXIMUM DIMENSIONS                                
      PARAMETER(Nmax=200,NFAC=20) 
      PARAMETER(Mmax=NFAC*Nmax) 
*                                                                              
      DIMENSION XL(Nmax),XR(Nmax),MF(Nmax),IFX(Nmax),XINIT(Nmax)
      DIMENSION S(Nmax*Mmax),F(Mmax) 
      DIMENSION X(Nmax),XMAX(Nmax),XMIN(Nmax),Y(Nmax) 
      DIMENSION IRN(0:Nmax)
      DIMENSION ZU(Nmax),WN(Nmax),CW(Nmax)
      LOGICAL CHECK                           

* Arrange the Panel Input
      call change(MERR)
      if (merr.ne.0) then
        call setcod(MERR)
        return
      endif

      call getpi('NOC',noc)
      call getpr('EPS',eps)
      call getpr('OMEGA',bign)
      call getpi('CANCEL',icanc)
      if(icanc.eq.0) then
        call setpi('CANCEL',1,IERR)
        return
      endif
* Import the dimension N (Number of Variables)
      call getdim(n,mmm)
      m= n*nfac
* Import the fix status.
      call getfix(IFX)
* Import current point
      call getx(XINIT)
* Set the lower and upper bounds for the variables.                            
* XL(i) i=1,...,N  are the lower bounds (Left margins)                         
* XR(i) i=1,...,N  are the upper bounds (Right margins)                        
      call getmar(mf,xl,xr)
*========================================================================
* STEP 0                                                                       
*========================================================================
* Initialize the Set S by sampling randomly M points.                          
      K = 0            
      NOE =0
      NOSE = 0
      CALL RSAMP(S,N,M,XL,XR)  
* Check if there are fixed variables. In that case reset the 
* variable to its specific initial value.
      DO 2 I = 1,N
         IF (IFX(I).eq.0) then
           DO 3 J=1,M
            S(I+N*J-N) = XINIT(I)
3          CONTINUE
         ENDIF
2     CONTINUE
* Evaluate the function at each point and collect the values.                  
      DO 4 I = 1,M            
	CALL ASSIGN(X,S,I,N,M) 
        F(I) = ACSQ(X,N)     
4     CONTINUE              
*                                                                             
*========================================================================
* STEP 1                                                                       
*========================================================================
*                                                                              
101   CONTINUE                                                  
      CALL  MINMAX(F,M,FMIN,IN,FMAX,IX)          
      IF (K .EQ. 0) THEN                        
        FMAX0 = FMAX                         
        FMIN0 = FMIN                          
      ENDIF                                 
      DO 10 I = 1,N                        
        XMAX(I) = S(I+N*IX-N)             
        XMIN(I) = S(I+N*IN-N)            
10    CONTINUE                          
      TEST = FMAX - FMIN               
      IF (TEST .LT. EPS .OR. NOE.GE.NOC) THEN        
       call setx(XMIN)
       call setval(fmin)
       return                                                       
      ENDIF                                                          
*                                                                              
*========================================================================
* Step 2                                                                       
*========================================================================
*                                                                              
* Pick N+1 integers from 1 to M randomly.                                      
* Store them in IRN(0:N)                                                       
*                                                                              
102   CONTINUE                                                  
      CALL  INTRAN(M,N,IRN)                                      
      PHI = BIGN*(FMAX-FMIN)**2/(FMAX0-FMIN0)                    
      DO 20 I=1,N                                              
       IP = IRN(I)                                             
       ZU(I) = 1./(F(IP)-FMIN+PHI)                                
20    CONTINUE                                    
      CALL NORMAL(ZU,N,WN)                                        
      DO 21 I = 1,N                                               
       CW(I) = 0                                               
        DO 22 J=1,N                                              
           CW(I) = CW(I) + WN(J)*S(I+N*IRN(J)-N) 
22      CONTINUE                                                 
21    CONTINUE                                                     
      FW = 0                                                       
      DO 23 I=1,N                                                   
       FW = FW + WN(I)*F(IRN(I))                                
23    CONTINUE                                                      
      I0 = IRN(0)                                                   
      FAC = (F(I0)-FW)/(FMAX-FMIN+PHI)                              
      IF (FW .LE. F(I0)) THEN                              
       DO 24 I=1,N                                                
         X(I) = (S(I+N*I0-N)-CW(I))*FAC+2*CW(I)-S(I+N*I0-N)       
24     CONTINUE                                                                
      ELSE                                                                     
       DO 25 I=1,N                                                             
         X(I) = (S(I+N*I0-N)-CW(I))*FAC+2*S(I+N*I0-N)-CW(I)                            
25     CONTINUE                                                                
      ENDIF                                                                    
      CALL DOMAIN(X,N,XL,XR,CHECK)                                             
      IF (CHECK) THEN                                                          
        FC = ACSQ(X,N)                                                       
        NOE = NOE + 1                                                           
	IF (NOE .GE. NOC) GO TO 101
      ELSE                                                                     
        GO TO 102                                                              
      ENDIF                                                                    
*                                                                              
*========================================================================
* STEP 3                                                                       
*========================================================================
*                                                                              
103   CONTINUE         
      IF ( FC .GE. FMAX ) THEN     
          RATE = DBLE(NOSE)/NOE   
          IF (RATE .GT. 0.5D0) THEN         
            K = K+1   
            GO TO 102       
          ENDIF      

          DO 30 I = 1,N         
            Y(I) = (CW(I)+S(I+N*IRN(N)-N) )/2.  
30        CONTINUE                    
          FY = ACSQ(Y,N)                                     
          NOE = NOE + 1          
          IF (FY .GE.FMAX) THEN                                    
            K = K+1                                                     
            GO TO 102                
          ELSE                                                          
            NOSE = NOSE + 1                                    
            CALL REPLAC(S,N,M,IX,Y)                                 
            F(IX)=FY                              
            K = K+1                 
            GO TO 101                                  
          ENDIF                 
      ENDIF  
      NOSE = NOSE + 1                       
      CALL REPLAC(S,N,M,IX,X)        
      F(IX)=FC              
      K = K+1              
      GO TO 101         
      END           

      SUBROUTINE ASSIGN(X,S,IP,N,M)                                            
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                                      
      DIMENSION S(N,M),X(N)                                                    
      DO 3 J=1,N                                                               
         X(J)=S(J,IP)                                                          
3     CONTINUE                                                                 
      END                                                                      
                                                                               
      SUBROUTINE DOMAIN(X,N,XL,XR,CHECK)                                       
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                                      
* Checks if point x is inside the domain.                                      
* I.e. checks if  x(i) is contained in [xl(i), xr(i)], i=1,2,..,n              
* INPUT:  x,n,xl,xr                                                            
* X       The input array to be checked                                        
* N       The elements of X to be considered. X(1) through X(N)                
* XL      The array containing the lower bounds. (Left margins)                
* XR      The array containing the upper bounds. (Right margins)               
* OUTPUT: CHECK                                                                
* CHECK   Logical variable. TRUE if X belongs to the domain. FALSE otherwise.  
      LOGICAL CHECK                                                            
      DIMENSION X(N),XL(N),XR(N)                                               
                                                                               
      CHECK = .TRUE.                                                           
      DO 1 I = 1,N                                                             
         IF (X(I) .LT. XL(I)) THEN                                             
	    CHECK = .FALSE.                                                           
	    RETURN                                                                    
	 ENDIF                                                                        
	 IF (X(I) .GT. XR(I)) THEN                                                    
	    CHECK = .FALSE.                                                           
	    RETURN                                                                    
         ENDIF                                                                 
 1    CONTINUE                                                                 
      END                                                                      
       SUBROUTINE MINMAX(F,M,FMIN,IN,FMAX,IX)                                  
       IMPLICIT DOUBLE PRECISION (A-H,O-Z)                                     
* Calculates the minimum and the maximum numbers                               
* contained in array f, and the corresponding indices.                         
* INPUT: F,M                                                                   
* F    The input array                                                         
* M    The number of elements to be considered. F(1) through F(M).             
* OUTPUT: FMIN,IN,FMAX,IX                                                      
* FMIN The minimum of the numbers in F                                         
* IN   The location of FMIN inside the array F. (FMIN = F(IN))                 
* FMAX The maximum of the numbers in F                                         
* IX   The location of FMAX inside the array F. (FMAX = F(IX))                 
* In case of multiple minima or maxima, the lowest indices                     
* are returned for IN and IX.                                                  
                                                                               
      DIMENSION F(*)                                                           
                                                                               
      IN = 1                                                                   
      IX = 1                                                                   
      FMIN = F(IN)                                                             
      FMAX = F(IX)    
                                                                               
      DO 1 I = 2,M                                                             
         FC = F(I)                                                             
	 IF (FC .LT. FMIN ) THEN                                                      
	     FMIN = FC                                                                
	     IN = I                                                                   
	 ENDIF                                                                        
	 IF (FC .GT. FMAX) THEN                                                       
	     FMAX = FC                                                                
	     IX = I                                                                   
	 ENDIF                                                                        
 1    CONTINUE                                                                 
      END                                                                      

      SUBROUTINE NORMAL(ZU,N,WN)                                               
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                                      
* Normalizes the array elements ZU(i)                                          
* Normalization here means: WN(i) = ZU(i)/{Sum_over_j=1,N ZU(j)}               
      DIMENSION ZU(*),WN(*)                                                    
                                                                               
      DEN = 0.                                                                 
      DO 1 I=1,N                                                               
        DEN = DEN + ZU(I)                                                      
1     CONTINUE                                                                 
                                                                               
      DO 2 I = 1,N                                                             
        WN(I) = ZU(I)/DEN                                                      
2     CONTINUE                                                                 
      END                                                                      
                                                                               
      SUBROUTINE INTRAN(M,N,IRN)                                               
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                                      
* Picks N+1 integers out of M at random.
* These integers are returned in IRN(0:N)
      PARAMETER (MMAX = 1000)                                                  
      DIMENSION IRN(0:N),INN(MMAX)                                             
      IF (N .LT. 1) THEN                                                       
       WRITE(*,*) 'IMPROPER VALUE:  N = ',N                                    
       STOP                                                                    
      ENDIF                                                                    
      IF (N+1 .GT. M) THEN                                                     
        WRITE(*,*) 'N CANNOT EXCEED: ', M-1                                    
	STOP                                                                          
      ENDIF                                                                    
      IF (M .GT. MMAX ) THEN                                                   
        WRITE(*,*) 'ROUTINE: INTRAN'                                           
        WRITE(*,*) 'INCREASE PARAMETER MMAX TO: ', M                           
	STOP                                                                          
      ENDIF                                                                    
      DO 1 I = 1,M                                                             
        INN(I) = I                                                             
 1    CONTINUE                                                                 
      IPOP = M                                                                 
      DO 2 I = 0,N                                                             
       XSI = RANM()                                                            
       IR = 1 + INT(IPOP*XSI)                                                  
       IRN(I) = INN(IR)                                                        
       IPOP = IPOP-1                                                           
       DO 20 J = IR,IPOP                                                       
       INN(J) = INN(J+1)                                                       
 20    CONTINUE                                                                
 2    CONTINUE                                                                 
      END                                                                      
                                                                               
C  ---------------------------------------------------------------------       
      SUBROUTINE REPLAC(S,N,M,JC,P)                                            
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                                      
* Replaces in the N x M matrix S, the JC column (N rows)                       
* by the N-vector P.                                                           
      DIMENSION S(N,M),P(N)                                                    
      DO 1 I = 1,N                                                             
        S(I, JC) = P(I)                                                        
1     CONTINUE                                                                 
      END                                                                      
                                                                               
      SUBROUTINE RSAMP(S,N,M,XL,XR)                                            
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                                      
* This routine samples at random M points in the N-dimensional                 
* space, inside the N-dimensional box defined by the bounds                    
* XL (lower) and XR (upper).                                                   
      DIMENSION S(N,M)                                                         
      DIMENSION XL(N),XR(N)                                                    
                                                                               
      DO 1 IP = 1,M                                                            
        DO 2 I = 1,N                                                           
	  S(I,IP)=RANDI(XL(I),XR(I))                                                  
 2      CONTINUE                                                               
 1    CONTINUE                                                                 
      END                                                                      

      SUBROUTINE WAVER(WN,N,Y,YAV)                                             
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                                      
      DIMENSION WN(N),Y(N)                                                     
                                                                               
      YAV = 0                                                                  
      DO 1 I=1,N                                                               
        YAV = YAV + WN(I)*Y(I)                                                 
1     CONTINUE                                                                 
      END                                                                      

