      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(*),XL(*),XR(*)                                               
                                                                               
      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                                                                      
      PROGRAM PRCGLO                                                           
* This is an implementation of the improved Price algorithm                    
* for global optimization.                                                     
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                                      
* INPUT PARAMETERS:                                                            
* N is the problem dimensionality
* M is the random sample size (Number of random points)
* EPS is the termination tolerance
* BIGN is a weigthing factor
* NOC is a maximum number of calls allowed.
      PARAMETER(N=10)                                                           
      PARAMETER(M=20*N)                                                        
      PARAMETER(EPS = 10.D-4)                                                  
      PARAMETER(BIGN= 1000.D0)                                                 
      PARAMETER(NOC = 50000)
*                                                                              
      DIMENSION XL(N),XR(N)                                                    
      DIMENSION S(N,M),F(M)                                                    
      DIMENSION X(N),XMAX(N),XMIN(N),Y(N)                                      
      DIMENSION IRN(0:N)                                                       
      DIMENSION ZU(N),WN(N),CW(N)                                              
      LOGICAL CHECK                                                            
                                                                               
* STEP 0                                                                       
* 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)                        
C One should prepare a file with filename: bounds.dat                          
C with the left and right margins in the following format:                     
C The file must contain N lines.                                               
C Each line contains two numbers, one for the lower and one for the            
C upper bound correspondingly.                                                 
      OPEN(1,FILE='bounds.dat')                                                
      DO 1 I = 1,N                                                             
       READ(1,*)XL(I),XR(I)                                                    
1     CONTINUE                                                                 
      CLOSE(1)                                                                 
* Initialize the Set S by sampling randomly M points.                          
      K = 0                                                                    
      NOE =0
      NOSE = 0
      CALL RSAMP(S,N,M,XL,XR)                                                  
* Evaluate the function at each point and collect the values.                  
      DO 2 I = 1,M                                                             
	CALL ASSIGN(X,S,I,N,M)                                                        
        F(I) = FUNMIN(X,N)                                                     
2     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,IX)                                                      
	XMIN(I) = S(I,IN)                                                             
10    CONTINUE                                                                 
      TEST = FMAX - FMIN                                                       
      OPEN(1,FILE='INIT')
      IF (TEST .LT. EPS .OR. NOE .ge. NOC) THEN            
       DO 11 I = 1,N                                                   
        WRITE(*,*)XMIN(I)                                                      
        WRITE(1,*)XMIN(I)                                                      
11     CONTINUE                                                                
       close(1)
       WRITE(*,*)'VALUE = ',FMIN                                               
       WRITE(*,*)'NUMBER OF EVALUATIONS: ',NOE+M                            
       STOP                                                                    
      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,IRN(J))                                           
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,I0)-CW(I))*FAC+2*CW(I)-S(I,I0)                            
24     CONTINUE                                                                
      ELSE                                                                     
       DO 25 I=1,N                                                             
         X(I) = (S(I,I0)-CW(I))*FAC+2*S(I,I0)-CW(I)                            
25     CONTINUE                                                                
      ENDIF                                                                    
      CALL DOMAIN(X,N,XL,XR,CHECK)                                             
      IF (CHECK) THEN                                                          
        FC = FUNMIN(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,IRN(N)) )/2.                                            
30        CONTINUE                                                             
          FY = FUNMIN(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 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)                                      
      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 RANLUX(RVEC,LENV)                                             
C  ---------------------------------------------------------------------       
C         Subtract-and-borrow random number generator proposed by              
C         Marsaglia and Zaman, implemented by F. James with the name           
C         RCARRY in 1991, and later improved by Martin Luescher                
C         in 1993 to produce "Luxury Pseudorandom Numbers".                    
C     Fortran 77 coded by F. James, 1993                                       
C                                                                              
C       references:                                                            
C  M. Luscher, Computer Physics Communications  79 (1994) 100                  
C  F. James, Computer Physics Communications 79 (1994) 111                     
C                                                                              
C   LUXURY LEVELS.                                                             
C   ------ ------      The available luxury levels are:                        
C                                                                              
C  level 0  (p=24): equivalent to the original RCARRY of Marsaglia             
C           and Zaman, very long period, but fails many tests.                 
C  level 1  (p=48): considerable improvement in quality over level 0,          
C           now passes the gap test, but still fails spectral test.            
C  level 2  (p=97): passes all known tests, but theoretically still            
C           defective.                                                         
C  level 3  (p=223): DEFAULT VALUE.  Any theoretically possible                
C           correlations have very small chance of being observed.             
C  level 4  (p=389): highest possible luxury, all 24 bits chaotic.             
C                                                                              
C!!! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++        
C!!!  Calling sequences for RANLUX:                                  ++        
C!!!      CALL RANLUX (RVEC, LEN)   returns a vector RVEC of LEN     ++        
C!!!                   32-bit random floating point numbers between  ++        
C!!!                   zero (not included) and one (also not incl.). ++        
C!!!      CALL RLUXGO(LUX,INT,K1,K2) initializes the generator from  ++        
C!!!               one 32-bit integer INT and sets Luxury Level LUX  ++        
C!!!               which is integer between zero and MAXLEV, or if   ++        
C!!!               LUX .GT. 24, it sets p=LUX directly.  K1 and K2   ++        
C!!!               should be set to zero unless restarting at a break++        
C!!!               point given by output of RLUXAT (see RLUXAT).     ++        
C!!!      CALL RLUXAT(LUX,INT,K1,K2) gets the values of four integers++        
C!!!               which can be used to restart the RANLUX generator ++        
C!!!               at the current point by calling RLUXGO.  K1 and K2++        
C!!!               specify how many numbers were generated since the ++        
C!!!               initialization with LUX and INT.  The restarting  ++        
C!!!               skips over  K1+K2*E9   numbers, so it can be long.++        
C!!!   A more efficient but less convenient way of restarting is by: ++        
C!!!      CALL RLUXIN(ISVEC)    restarts the generator from vector   ++        
C!!!                   ISVEC of 25 32-bit integers (see RLUXUT)      ++        
C!!!      CALL RLUXUT(ISVEC)    outputs the current values of the 25 ++        
C!!!                 32-bit integer seeds, to be used for restarting ++        
C!!!      ISVEC must be dimensioned 25 in the calling program        ++        
C!!! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++        
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                                      
      DIMENSION RVEC(LENV)                                                     
      DIMENSION SEEDS(24), ISEEDS(24), ISDEXT(25)                              
      PARAMETER (MAXLEV=4, LXDFLT=3)                                           
      DIMENSION NDSKIP(0:MAXLEV)                                               
      DIMENSION NEXT(24)                                                       
      PARAMETER (TWOP12=4096., IGIGA=1000000000,JSDFLT=314159265)              
      PARAMETER (ITWO24=2**24, ICONS=2147483563)                               
      SAVE NOTYET, I24, J24, CARRY, SEEDS, TWOM24, TWOM12, LUXLEV              
      SAVE NSKIP, NDSKIP, IN24, NEXT, KOUNT, MKOUNT, INSEED                    
      INTEGER LUXLEV                                                           
      LOGICAL NOTYET                                                           
      DATA NOTYET, LUXLEV, IN24, KOUNT, MKOUNT /.TRUE., LXDFLT, 0,0,0/         
      DATA I24,J24,CARRY/24,10,0./                                             
C                               default                                        
C  Luxury Level   0     1     2   *3*    4                                     
      DATA NDSKIP/0,   24,   73,  199,  365 /                                  
Corresponds to p=24    48    97   223   389                                    
C     time factor 1     2     3     6    10   on slow workstation              
C                 1    1.5    2     3     5   on fast mainframe                
C                                                                              
C  NOTYET is .TRUE. if no initialization has been performed yet.               
C              Default Initialization by Multiplicative Congruential           
      IF (NOTYET) THEN                                                         
         NOTYET = .FALSE.                                                      
         JSEED = JSDFLT                                                        
         INSEED = JSEED                                                        
CC         WRITE(6,'(A,I12)') ' RANLUX DEFAULT INITIALIZATION: ',JSEED         
         LUXLEV = LXDFLT                                                       
         NSKIP = NDSKIP(LUXLEV)                                                
         IN24 = 0                                                              
         KOUNT = 0                                                             
         MKOUNT = 0                                                            
CC         LP = NSKIP + 24                                                     
CC         WRITE(6,'(A,I2,A,I4)')  ' RANLUX DEFAULT LUXURY LEVEL =  ',         
CC     +        LUXLEV,'      p =',LP                                          
            TWOM24 = 1.                                                        
         DO 25 I= 1, 24                                                        
            TWOM24 = TWOM24 * 0.5                                              
         K = JSEED/53668                                                       
         JSEED = 40014*(JSEED-K*53668) -K*12211                                
         IF (JSEED .LT. 0)  JSEED = JSEED+ICONS                                
         ISEEDS(I) = MOD(JSEED,ITWO24)                                         
   25    CONTINUE                                                              
         TWOM12 = TWOM24 * 4096.                                               
         DO 50 I= 1,24                                                         
         SEEDS(I) = REAL(ISEEDS(I))*TWOM24                                     
         NEXT(I) = I-1                                                         
   50    CONTINUE                                                              
         NEXT(1) = 24                                                          
         I24 = 24                                                              
         J24 = 10                                                              
         CARRY = 0.                                                            
         IF (SEEDS(24) .EQ. 0.) CARRY = TWOM24                                 
      ENDIF                                                                    
C                                                                              
C          The Generator proper: "Subtract-with-borrow",                       
C          as proposed by Marsaglia and Zaman,                                 
C          Florida State University, March, 1989                               
C                                                                              
      DO 100 IVEC= 1, LENV                                                     
      UNI = SEEDS(J24) - SEEDS(I24) - CARRY                                    
      IF (UNI .LT. 0.)  THEN                                                   
         UNI = UNI + 1.0                                                       
         CARRY = TWOM24                                                        
      ELSE                                                                     
         CARRY = 0.                                                            
      ENDIF                                                                    
      SEEDS(I24) = UNI                                                         
      I24 = NEXT(I24)                                                          
      J24 = NEXT(J24)                                                          
      RVEC(IVEC) = UNI                                                         
C  small numbers (with less than 12 "significant" bits) are "padded".          
      IF (UNI .LT. TWOM12)  THEN                                               
         RVEC(IVEC) = RVEC(IVEC) + TWOM24*SEEDS(J24)                           
C        and zero is forbidden in case someone takes a logarithm               
         IF (RVEC(IVEC) .EQ. 0.)  RVEC(IVEC) = TWOM24*TWOM24                   
      ENDIF                                                                    
C        Skipping to luxury.  As proposed by Martin Luscher.                   
      IN24 = IN24 + 1                                                          
      IF (IN24 .EQ. 24)  THEN                                                  
         IN24 = 0                                                              
         KOUNT = KOUNT + NSKIP                                                 
         DO 90 ISK= 1, NSKIP                                                   
         UNI = SEEDS(J24) - SEEDS(I24) - CARRY                                 
         IF (UNI .LT. 0.)  THEN                                                
            UNI = UNI + 1.0                                                    
            CARRY = TWOM24                                                     
         ELSE                                                                  
            CARRY = 0.                                                         
         ENDIF                                                                 
         SEEDS(I24) = UNI                                                      
         I24 = NEXT(I24)                                                       
         J24 = NEXT(J24)                                                       
   90    CONTINUE                                                              
      ENDIF                                                                    
  100 CONTINUE                                                                 
      KOUNT = KOUNT + LENV                                                     
      IF (KOUNT .GE. IGIGA)  THEN                                              
         MKOUNT = MKOUNT + 1                                                   
         KOUNT = KOUNT - IGIGA                                                 
      ENDIF                                                                    
      RETURN                                                                   
C                                                                              
C           Entry to input and float integer seeds from previous run           
      ENTRY RLUXIN(ISDEXT)                                                     
*     IF block added by Phillip Helbig after correpondence with James          
      IF (NOTYET) THEN                                                         
CC        WRITE(6,'(A)')  ' PROPER RESULTS ONLY WITH INITIALISATION FROM       
CC    $25 INTEGERS OBTAINED WITH RLUXUT'                                       
         NOTYET = .FALSE.                                                      
      ENDIF                                                                    
         TWOM24 = 1.                                                           
         DO 195 I= 1, 24                                                       
         NEXT(I) = I-1                                                         
  195    TWOM24 = TWOM24 * 0.5                                                 
         NEXT(1) = 24                                                          
         TWOM12 = TWOM24 * 4096.                                               
CC    WRITE(6,'(A)') ' FULL INITIALIZATION OF RANLUX WITH 25 INTEGERS:'        
CC    WRITE(6,'(5X,5I12)') ISDEXT                                              
      DO 200 I= 1, 24                                                          
      SEEDS(I) = REAL(ISDEXT(I))*TWOM24                                        
  200 CONTINUE                                                                 
      CARRY = 0.                                                               
      IF (ISDEXT(25) .LT. 0)  CARRY = TWOM24                                   
      ISD = IABS(ISDEXT(25))                                                   
      I24 = MOD(ISD,100)                                                       
      ISD = ISD/100                                                            
      J24 = MOD(ISD,100)                                                       
      ISD = ISD/100                                                            
      IN24 = MOD(ISD,100)                                                      
      ISD = ISD/100                                                            
      LUXLEV = ISD                                                             
        IF (LUXLEV .LE. MAXLEV) THEN                                           
          NSKIP = NDSKIP(LUXLEV)                                               
CC        WRITE (6,'(A,I2)') ' RANLUX LUXURY LEVEL SET BY RLUXIN TO: ',        
CC   +                         LUXLEV                                          
        ELSE  IF (LUXLEV .GE. 24) THEN                                         
          NSKIP = LUXLEV - 24                                                  
CC        WRITE (6,'(A,I5)') ' RANLUX P-VALUE SET BY RLUXIN TO:',LUXLEV        
        ELSE                                                                   
          NSKIP = NDSKIP(MAXLEV)                                               
CC        WRITE (6,'(A,I5)') ' RANLUX ILLEGAL LUXURY RLUXIN: ',LUXLEV          
          LUXLEV = MAXLEV                                                      
        ENDIF                                                                  
      INSEED = -1                                                              
      RETURN                                                                   
C                                                                              
C                    Entry to output seeds as integers                         
      ENTRY RLUXUT(ISDEXT)                                                     
      DO 300 I= 1, 24                                                          
         ISDEXT(I) = INT(SEEDS(I)*TWOP12*TWOP12)                               
  300 CONTINUE                                                                 
      ISDEXT(25) = I24 + 100*J24 + 10000*IN24 + 1000000*LUXLEV                 
      IF (CARRY .GT. 0.)  ISDEXT(25) = -ISDEXT(25)                             
      RETURN                                                                   
C                                                                              
C                    Entry to output the "convenient" restart point            
      ENTRY RLUXAT(LOUT,INOUT,K1,K2)                                           
      LOUT = LUXLEV                                                            
      INOUT = INSEED                                                           
      K1 = KOUNT                                                               
      K2 = MKOUNT                                                              
      RETURN                                                                   
C                                                                              
C                    Entry to initialize from one or three integers            
      ENTRY RLUXGO(LUX,INS,K1,K2)                                              
         IF (LUX .LT. 0) THEN                                                  
            LUXLEV = LXDFLT                                                    
         ELSE IF (LUX .LE. MAXLEV) THEN                                        
            LUXLEV = LUX                                                       
         ELSE IF (LUX .LT. 24 .OR. LUX .GT. 2000) THEN                         
            LUXLEV = MAXLEV                                                    
CC          WRITE (6,'(A,I7)') ' RANLUX ILLEGAL LUXURY RLUXGO: ',LUX           
         ELSE                                                                  
            LUXLEV = LUX                                                       
            DO 310 ILX= 0, MAXLEV                                              
              IF (LUX .EQ. NDSKIP(ILX)+24)  LUXLEV = ILX                       
  310       CONTINUE                                                           
         ENDIF                                                                 
      IF (LUXLEV .LE. MAXLEV)  THEN                                            
         NSKIP = NDSKIP(LUXLEV)                                                
CC       WRITE(6,'(A,I2,A,I4)') ' RANLUX LUXURY LEVEL SET BY RLUXGO :',        
CC   +        LUXLEV,'     P=', NSKIP+24                                       
      ELSE                                                                     
          NSKIP = LUXLEV - 24                                                  
CC        WRITE (6,'(A,I5)') ' RANLUX P-VALUE SET BY RLUXGO TO:',LUXLEV        
      ENDIF                                                                    
      IN24 = 0                                                                 
CC    IF (INS .LT. 0)  WRITE (6,'(A)')                                         
CC   +   ' Illegal initialization by RLUXGO, negative input seed'              
      IF (INS .GT. 0)  THEN                                                    
        JSEED = INS                                                            
CC      WRITE(6,'(A,3I12)') ' RANLUX INITIALIZED BY RLUXGO FROM SEEDS',        
CC   +      JSEED, K1,K2                                                       
      ELSE                                                                     
        JSEED = JSDFLT                                                         
CC      WRITE(6,'(A)')' RANLUX INITIALIZED BY RLUXGO FROM DEFAULT SEED'        
      ENDIF                                                                    
      INSEED = JSEED                                                           
      NOTYET = .FALSE.                                                         
      TWOM24 = 1.                                                              
         DO 325 I= 1, 24                                                       
           TWOM24 = TWOM24 * 0.5                                               
         K = JSEED/53668                                                       
         JSEED = 40014*(JSEED-K*53668) -K*12211                                
         IF (JSEED .LT. 0)  JSEED = JSEED+ICONS                                
         ISEEDS(I) = MOD(JSEED,ITWO24)                                         
  325    CONTINUE                                                              
      TWOM12 = TWOM24 * 4096.                                                  
         DO 350 I= 1,24                                                        
         SEEDS(I) = REAL(ISEEDS(I))*TWOM24                                     
         NEXT(I) = I-1                                                         
  350    CONTINUE                                                              
      NEXT(1) = 24                                                             
      I24 = 24                                                                 
      J24 = 10                                                                 
      CARRY = 0.                                                               
      IF (SEEDS(24) .EQ. 0.) CARRY = TWOM24                                    
C        If restarting at a break point, skip K1 + IGIGA*K2                    
C        Note that this is the number of numbers delivered to                  
C        the user PLUS the number skipped (if luxury .GT. 0).                  
      KOUNT = K1                                                               
      MKOUNT = K2                                                              
      IF (K1+K2 .NE. 0)  THEN                                                  
        DO 500 IOUTER= 1, K2+1                                                 
          INNER = IGIGA                                                        
          IF (IOUTER .EQ. K2+1)  INNER = K1                                    
          DO 450 ISK= 1, INNER                                                 
            UNI = SEEDS(J24) - SEEDS(I24) - CARRY                              
            IF (UNI .LT. 0.)  THEN                                             
               UNI = UNI + 1.0                                                 
               CARRY = TWOM24                                                  
            ELSE                                                               
               CARRY = 0.                                                      
            ENDIF                                                              
            SEEDS(I24) = UNI                                                   
            I24 = NEXT(I24)                                                    
            J24 = NEXT(J24)                                                    
  450     CONTINUE                                                             
  500   CONTINUE                                                               
C         Get the right value of IN24 by direct calculation                    
        IN24 = MOD(KOUNT, NSKIP+24)                                            
        IF (MKOUNT .GT. 0)  THEN                                               
           IZIP = MOD(IGIGA, NSKIP+24)                                         
           IZIP2 = MKOUNT*IZIP + IN24                                          
           IN24 = MOD(IZIP2, NSKIP+24)                                         
        ENDIF                                                                  
C       Now IN24 had better be between zero and 23 inclusive                   
        IF (IN24 .GT. 23) THEN                                                 
CC         WRITE (6,'(A/A,3I11,A,I5)')                                         
CC   +    '  Error in RESTARTING with RLUXGO:','  The values', INS,            
CC   +     K1, K2, ' cannot occur at luxury level', LUXLEV                     
           IN24 = 0                                                            
        ENDIF                                                                  
      ENDIF                                                                    
      RETURN                                                                   
      END                                                                      
C  ---------------------------------------------------------------------       
      FUNCTION RANM ( )                                                        
C  ---------------------------------------------------------------------       
C                                                                              
C  Description:                                                                
C    Function RANM returns a random number, uniformly distributed              
C    in the interval (0,1). It uses the RANLUX random number generator         
C    due to  Luscher and James. For details see:                               
C      Luscher M., Comput. Phys. Commun. 79 (1994) 100.                        
C      James F., Comput. Phys. Commun. 79 (1994) 111.                          
C                                                                              
C  ---------------------------------------------------------------------       
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                                      
C  Local variables:                                                            
      PARAMETER ( LVEC = 1 )                                                   
      DIMENSION RVEC(LVEC)                                                     
C  ---------------------------------------------------------------------       
C                                                                              
      CALL RANLUX(RVEC,LVEC)                                                   
      RANM = RVEC(1)                                                           
C                                                                              
      END                                                                      
C  ---------------------------------------------------------------------       
      FUNCTION RANDI ( A, B )                                                  
C  ---------------------------------------------------------------------       
C                                                                              
C  Description:                                                                
C   Function RANDI returns a random number between A and B.                    
C                                                                              
C  ---------------------------------------------------------------------       
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                                      
C                                                                              
      R = RANM()                                                               
      RANDI = A+(B-A)*R                                                        
C                                                                              
      END                                                                      
      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,*)                                                         
      DIMENSION XL(*),XR(*)                                                    
                                                                               
      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(*),Y(*)                                                     
                                                                               
      YAV = 0                                                                  
      DO 1 I=1,N                                                               
        YAV = YAV + WN(I)*Y(I)                                                 
1     CONTINUE                                                                 
      END                                                                      

      function funmin(x,n)
      implicit double precision (a-h,o-z)
      parameter (pi=3.14159265359d0)
      dimension x(n)
      funmin = 0
      do 1 i=1,n-1
      funmin = funmin + (x(i)-1)**2*(1+10*sin(pi*x(i+1))**2)
1     continue
      funmin = funmin + 10*sin(pi*x(1))**2+(x(n)-1)**2
      funmin = funmin*pi/n
      end
