C ======================================================================
      PROGRAM GENOMINO
C ======================================================================
      implicit double precision (a-h,o-z)
      parameter (n= 2, m= 200, ng =  800)
      parameter (delpc =   5.d-2, selpc=  5.d-2,mutpc=  5.d-2)
      parameter (PMUTE=  5.d0)
      dimension POP(n,m),S(n,m),X(n),FVAL(m),PROB(m),WA(n)      
      dimension XL(n),XR(n),PA1(n),PA2(n),CH1(n),CH2(n)
      dimension INDX(m),SCR(m)
      data XL/(n)*-5.d0/
      data XR/(n)*5.d0/
      ndel = max(1.d0,delpc*m)
      nsel = max(1.d0,selpc*m)
      nmut = max(1.d0,mutpc*m)
      ncro = m-ndel-nsel-nmut
C  NCRO must be even. If not decrement it by one and increment NMUT
      if ((ncro/2)*2 .ne. ncro) Then
              ncro = ncro -1
              nmut = nmut +1
      endif        
      write(*,'(57(''=''))')
      write(*,*)'GENOMINO RUNS WITH THE FOLLOWING DATA:'
      write(*,*)'Number of Variables: ', n
      write(*,*) 'D-Elitism number: ',ndel
      write(*,*) 'P-Elitism number: ',nsel
      write(*,*) 'Crossover number: ',ncro
      write(*,*) 'Mutation  number: ',nmut
      write(*,*) 'Mutation    Rate: ',pmute
      write(*,*) 'Population  Size: ',m
      write(*,*) 'Max. Generations: ',ng 
      write(*,'(57(''=''))')
C  Create initial generation
      CALL CREATE(M,POP,N,XL,XR,FVAL,WA)
      IGENER = 1
30    CONTINUE
      if (IGENER.eq.ng) go to 31 
      IGENER = IGENER + 1
*  Calculate probabilities from  FVAL
      CALL calpro(fval,m,PROB) 
*  Find the positions of the NDEL lowest values in FVAL, and store them in
*   array INDX(M). 
      CALL ELITE(FVAL,M,INDX,NDEL,SCR)
      do 21 i =1,ndel
         iel = indx(i)
         CALL PICKG(POP,N,M,IEL,WA)
*  Keep the best NDEL points in the next generation.
*  (Deterministic  Elitism)
*  S(n,m) is the next generation population.
         CALL SETCOL(S,N,M,I,WA)
21    continue         
      IFIL = ndel
*  Selection Operation. (Roulette Wheel)
*  Select for the next generation NSEL points, probabilistically.
*  Points corresponding to lower values of the objective, have higher 
*  probability to "survive". (Probabilistic Elitism).
      do 1 is = 1,nsel
         IFIL = IFIL + 1
         CALL select(prob,m,j) 
         CALL PICKG(POP,N,M,J,WA)
*  Start forming the next generation population S.
         CALL SETCOL(S,N,M,IFIL,WA)
1     continue
      do 2 ic = 1,ncro/2
         CALL select(PROB,M,J1)
*  To avoid identical parents, temporarily set PROB(J1)=0.
         temp = PROB(j1)
         PROB(J1)=0.
         CALL select(PROB,M,J2)
*  Reset PROB(J1) back to its true value.
         PROB(J1) = temp
*  Construct the parents PA1, PA2
         CALL PICKG(POP,N,M,J1,PA1)
         CALL PICKG(POP,N,M,J2,PA2)
*  Crossover Operation: Create the children CH1, CH2
*             CALL crossd(pa1,pa2,n,ch1,ch2) 
              CALL crossi(pa1,pa2,n,ch1,ch2,xl,xr) 
*  Continue filling the next generation population S.
         IFIL = IFIL + 1
         CALL SETCOL(S,N,M,IFIL,CH1)
         IFIL = IFIL + 1
         CALL SETCOL(S,N,M,IFIL,CH2)
2     CONTINUE
*  Mutation Operation
         do 3 j=1,nmut
           IFIL = IFIL + 1
           CALL select(PROB,M,J1)
           CALL PICKG(POP,N,M,J1,PA1)
           CALL MUTATE(XL,XR,PA1,n,IGENER,M,PMUTE,CH1)
           CALL SETCOL(S,N,M,IFIL,CH1)
3        continue
c  Start a new generation
      CALL equate(POP,N,M,S)
      CALL evaluate(POP,n,m,FVAL,wa)
            CALL MINEL(FVAL,M,IEL)
*           write(*,'(1x,5(g14.7,1x))') FVAL(iel),(S(j,iel),j=1,n) 
      go to 30
31    continue
      write(*,88)'GA:','FUNCTION CALLS','MINIMUM','POINT'
88    format(1x,a5,T9,a15,T32,a7,T40,A17)
      write(*,89) m*ng,FVAL(IEL),(S(j,iel),j=1,n)
            write(*,'(57(''-''))')
89    format(T13,i7,T29,D14.7,T45,d14.7,10(/t45,d14.7))
      end
C ======================================================================
      SUBROUTINE CREATE(M,S,N,XL,XR,FVAL,WA)
C ======================================================================
*  Creates a population of M members.
*  Each member is an array with N components (i.e. a point in R^N)
*  XL(I), XR(I) are the lower and upper bounds for the I-th component. 
*  S(N,M) contains the population, i.e M,  N-vectors.
*  FVAL(J) contains the value of the objective function evaluated at the 
*  point corresponding to J-th column of S.
*  WA is a work array. 
** INPUT:  M,N,XL,XR
** OUTPUT: S,FVAL
      implicit double precision (a-h,o-z)
      dimension s(n,m),xl(n),xr(n),wa(n),fval(m)
      do 1 i=1,m
         do 2 j=1,n
           s(j,i) = (xr(j)-xl(j))*ranm()+xl(j)
           wa(j) = s(j,i)
2        continue
           fval(i) = funmin(wa,n)
1     continue
      end
C ======================================================================
      subroutine crossd(p1,p2,n,c1,c2) 
C ======================================================================
*  Implements the discrete crossover operation, between two parents.
*  P1, P2 are the parent members.
*  C1, C2 are the children  members.
*      
      implicit double precision (a-h,o-z)
      dimension p1(n),p2(n),c1(n),c2(n)
      logical l1
      do 1 i=1,n
         xi = ranm()
         l1 = xi .le. 5.d-1 
         if (l1) then
                 c1(i) = p1(i)
                 c2(i) = p2(i)
         else
                 c1(i) = p2(i)
                 c2(i) = p1(i)
         endif
1     continue
      end         
C ======================================================================
      subroutine crossi(p1,p2,n,c1,c2,xl,xr) 
C ======================================================================
*  Implements the intermediate crossover operation, between two parents.
*  P1, P2 are the parent members.
*  C1, C2 are the children  members.
*  Care is taken that the children are inside the box.
      implicit double precision (a-h,o-z)
      parameter (d=0.25d0)
      dimension p1(*),p2(*),c1(*),c2(*),xl(*),xr(*)
      da = 1+2*d
      a0 = -d
      do 1 i=1,n
         alfa1 = da*ranm()+a0
         alfa2 = da*ranm()+a0
         c1(i) = p1(i) + alfa1*(p2(i)-p1(i))
           if (c1(i) .lt. xl(i) ) c1(i) = xl(i)
           if (c1(i) .gt. xr(i) ) c1(i) = xr(i)
         c2(i) = p1(i) + alfa2*(p2(i)-p1(i))
           if (c2(i) .lt. xl(i) ) c2(i) = xl(i)
           if (c2(i) .gt. xr(i) ) c2(i) = xr(i)
1     continue
      end
C ======================================================================
      subroutine calpro(fval,m,pr) 
C ======================================================================
*  Calculates the selection probabilities.
** INPUT:   FVAL, M
** OUTPUT:  PR
      implicit double precision (a-h,o-z)
      parameter (eps = 1.d-7)
      dimension fval(m),pr(m)
      sum = 0
      fmax = fval(1)
         do 1 i = 1,m
            sum = sum + fval(i)
            if (fval(i) .gt. fmax) fmax = fval(i)
1        continue
         div = m*fmax - sum
         if (div .eq. 0.d0) STOP 'CALPRO'
         do 2 i =1,m
            pr(i) = (fmax - fval(i))/div
            pr(i) = (pr(i)+eps)/(1+m*eps)
2        continue
      end
C ======================================================================
      subroutine select(pr,m,j) 
C ======================================================================
*  Implements the selection operation.
*  Selects according to the probabilities ( PR(M) ) one point.
**  INPUT: PR, M
**  OUTPUT: J, the index of the selected point.

      implicit double precision (a-h,o-z)
      dimension pr(m)
      xi = ranm()
      seg = 0.d0
          do 1 i = 1,m
             seg = seg + pr(i)
             if (xi .le. seg) Then
                     j=i
                     return
             endif
1         continue
      end
C ======================================================================
      subroutine equate(POP,N,M,S)
C ======================================================================
      implicit double precision (a-h,o-z)
      dimension POP(n,m), S(n,m)
      do 1 i = 1,n
         do 2 j = 1,m
             POP(i,j) = S(i,j)
2        continue
1     continue
      end
C ======================================================================
      subroutine evaluate(POP,n,m,FVAL,wa)
C ======================================================================
      implicit double precision (a-h,o-z)
      dimension POP(n,m),FVAL(m),wa(n)
      do 1 i=1,m
         do 2 j=1,n
           wa(j) = POP(j,i)
2        continue
           fval(i) = funmin(wa,n)
1     continue
      end
C ======================================================================
      subroutine PICKG(POP,N,M,ICL,P)
C ======================================================================
*  Selects the ICL column of POP, and stores it to P
      implicit double precision (a-h,o-z)
      dimension POP(n,m),P(n) 
      do 1 i=1,n
         P(i) = POP(i,ICL)
1     continue
      end
C ======================================================================
      subroutine SETCOL(POP,N,M,ICL,P)
C ======================================================================
*  Assigns P, to the ICL column of POP.
      implicit double precision (a-h,o-z)
      dimension POP(n,m),P(n) 
      do 1 i=1,n
         POP(i,ICL)=P(i)
1     continue
      end
C ======================================================================
      subroutine MUTATE(XL,XR,PA1,n,IGENER,M,PMUTE,CH1)
C ======================================================================
      implicit double precision (a-h,o-z)
      dimension PA1(n),CH1(n),XL(n),XR(n)
      b = (1.d0-dble(IGENER)/dble(M))*PMUTE
      do 1 i=1,n
         xi = ranm()
         yi = ranm()
         if (xi .gt. 0.5d0) Then
            ch1(i) = pa1(i) + (xr(i)-pa1(i))*(1.d0-yi**b)
         else
            ch1(i) = pa1(i) + (pa1(i)-xl(i))*(1.d0-yi**b)
         endif
1     continue
      end
C ======================================================================
      subroutine ELITE(FVAL,M,INDX,K,SCR)
C ======================================================================
*  INPUT : M, FVAL(M),K
*  OUTPUT: INDX(M), the positions of the K smallest values in FVAL(M)
*  SCR(M) is a workspace array
      implicit double precision (a-h,o-z)
      dimension FVAL(M),INDX(M),SCR(M)
      fmax = fval(1)
      do 10 i=1,m
        scr(i) = fval(i)
        indx(i) = 0
        if (fmax .lt. fval(i))fmax = fval(i)
10    continue
      do 1 i=1,k
        call minel(scr,m,iel)
          indx(i) = iel
          scr(iel)= fmax
1     continue
      end
C ======================================================================
      subroutine MINEL(FVAL,M,IEL)
C ======================================================================
*  INPUT : M, FVAL(M)
*  OUTPUT: IEL, the position of the smallest value in FVAL(M)
      implicit double precision (a-h,o-z)
      dimension FVAL(M)
      fmin = fval(1)
      iel = 1
        do 1 i = 1,m
          if (fval(i) .le. fmin) Then
            fmin = fval(i)
            iel = i
          endif
1       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
C  ---------------------------------------------------------------------
