MEVSIM/ranlux.F

   1 C*
   2 C* $Id$
   3 C*
   4 C#include "gen/pilot.h"
   5       SUBROUTINE RANLUX(RVEC,LENV)
   6 C         Subtract-and-borrow random number generator proposed by
   7 C         Marsaglia and Zaman, implemented by F. James with the name
   8 C         RCARRY in 1991, and later improved by Martin Luescher
   9 C         in 1993 to produce "Luxury Pseudorandom Numbers".
  10 C     Fortran 77 coded by F. James, 1993
  11 C
  12 C   LUXURY LEVELS.
  13 C   ------ ------      The available luxury levels are:
  14 C
  15 C  level 0  (p=24): equivalent to the original RCARRY of Marsaglia
  16 C           and Zaman, very long period, but fails many tests.
  17 C  level 1  (p=48): considerable improvement in quality over level 0,
  18 C           now passes the gap test, but still fails spectral test.
  19 C  level 2  (p=97): passes all known tests, but theoretically still
  20 C           defective.
  21 C  level 3  (p=223): DEFAULT VALUE.  Any theoretically possible
  22 C           correlations have very small chance of being observed.
  23 C  level 4  (p=389): highest possible luxury, all 24 bits chaotic.
  24 C
  25 C!!! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  26 C!!!  Calling sequences for RANLUX:                                  ++
  27 C!!!      CALL RANLUX (RVEC, LEN)   returns a vector RVEC of LEN     ++
  28 C!!!                   32-bit random floating point numbers between  ++
  29 C!!!                   zero (not included) and one (also not incl.). ++
  30 C!!!      CALL RLUXGO(LUX,INT,K1,K2) initializes the generator from  ++
  31 C!!!               one 32-bit integer INT and sets Luxury Level LUX  ++
  32 C!!!               which is integer between zero and MAXLEV, or if   ++
  33 C!!!               LUX .GT. 24, it sets p=LUX directly.  K1 and K2   ++
  34 C!!!               should be set to zero unless restarting at a break++
  35 C!!!               point given by output of RLUXAT (see RLUXAT).     ++
  36 C!!!      CALL RLUXAT(LUX,INT,K1,K2) gets the values of four integers++
  37 C!!!               which can be used to restart the RANLUX generator ++
  38 C!!!               at the current point by calling RLUXGO.  K1 and K2++
  39 C!!!               specify how many numbers were generated since the ++
  40 C!!!               initialization with LUX and INT.  The restarting  ++
  41 C!!!               skips over  K1+K2*E9   numbers, so it can be long.++
  42 C!!!   A more efficient but less convenient way of restarting is by: ++
  43 C!!!      CALL RLUXIN(ISVEC)    restarts the generator from vector   ++
  44 C!!!                   ISVEC of 25 32-bit integers (see RLUXUT)      ++
  45 C!!!      CALL RLUXUT(ISVEC)    outputs the current values of the 25 ++
  46 C!!!                 32-bit integer seeds, to be used for restarting ++
  47 C!!!      ISVEC must be dimensioned 25 in the calling program        ++
  48 C!!! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  49       DIMENSION RVEC(LENV)
  50       DIMENSION SEEDS(24), ISEEDS(24), ISDEXT(25)
  51       PARAMETER (MAXLEV=4, LXDFLT=3)
  52       DIMENSION NDSKIP(0:MAXLEV)
  53       DIMENSION NEXT(24)
  54       PARAMETER (TWOP12=4096., IGIGA=1000000000,JSDFLT=314159265)
  55       PARAMETER (ITWO24=2**24, ICONS=2147483563)
  56       SAVE NOTYET, I24, J24, CARRY, SEEDS, TWOM24, TWOM12, LUXLEV
  57       SAVE NSKIP, NDSKIP, IN24, NEXT, KOUNT, MKOUNT, INSEED
  58       INTEGER LUXLEV
  59       LOGICAL NOTYET
  60       DATA NOTYET, LUXLEV, IN24, KOUNT, MKOUNT /.TRUE., LXDFLT, 0,0,0/
  61       DATA I24,J24,CARRY/24,10,0./
  62 C                               default
  63 C  Luxury Level   0     1     2   *3*    4
  64       DATA NDSKIP/0,   24,   73,  199,  365 /
  65 Corresponds to p=24    48    97   223   389
  66 C     time factor 1     2     3     6    10   on slow workstation
  67 C                 1    1.5    2     3     5   on fast mainframe
  68 C
  69 C  NOTYET is .TRUE. if no initialization has been performed yet.
  70 C              Default Initialization by Multiplicative Congruential
  71       IF (NOTYET) THEN
  72          NOTYET = .FALSE.
  73          JSEED = JSDFLT
  74          INSEED = JSEED
  75          WRITE(6,'(A,I12)') ' RANLUX DEFAULT INITIALIZATION: ',JSEED
  76          LUXLEV = LXDFLT
  77          NSKIP = NDSKIP(LUXLEV)
  78          LP = NSKIP + 24
  79          IN24 = 0
  80          KOUNT = 0
  81          MKOUNT = 0
  82          WRITE(6,'(A,I2,A,I4)')  ' RANLUX DEFAULT LUXURY LEVEL =  ',
  83      +        LUXLEV,'      p =',LP
  84             TWOM24 = 1.
  85          DO 25 I= 1, 24
  86             TWOM24 = TWOM24 * 0.5
  87          K = JSEED/53668
  88          JSEED = 40014*(JSEED-K*53668) -K*12211
  89          IF (JSEED .LT. 0)  JSEED = JSEED+ICONS
  90          ISEEDS(I) = MOD(JSEED,ITWO24)
  91    25    CONTINUE
  92          TWOM12 = TWOM24 * 4096.
  93          DO 50 I= 1,24
  94          SEEDS(I) = REAL(ISEEDS(I))*TWOM24
  95          NEXT(I) = I-1
  96    50    CONTINUE
  97          NEXT(1) = 24
  98          I24 = 24
  99          J24 = 10
 100          CARRY = 0.
 101          IF (SEEDS(24) .EQ. 0.) CARRY = TWOM24
 102       ENDIF
 103 C
 104 C          The Generator proper: "Subtract-with-borrow",
 105 C          as proposed by Marsaglia and Zaman,
 106 C          Florida State University, March, 1989
 107 C
 108       DO 100 IVEC= 1, LENV
 109       UNI = SEEDS(J24) - SEEDS(I24) - CARRY
 110       IF (UNI .LT. 0.)  THEN
 111          UNI = UNI + 1.0
 112          CARRY = TWOM24
 113       ELSE
 114          CARRY = 0.
 115       ENDIF
 116       SEEDS(I24) = UNI
 117       I24 = NEXT(I24)
 118       J24 = NEXT(J24)
 119       RVEC(IVEC) = UNI
 120 C  small numbers (with less than 12 "significant" bits) are "padded".
 121       IF (UNI .LT. TWOM12)  THEN
 122          RVEC(IVEC) = RVEC(IVEC) + TWOM24*SEEDS(J24)
 123 C        and zero is forbidden in case someone takes a logarithm
 124          IF (RVEC(IVEC) .EQ. 0.)  RVEC(IVEC) = TWOM24*TWOM24
 125       ENDIF
 126 C        Skipping to luxury.  As proposed by Martin Luscher.
 127       IN24 = IN24 + 1
 128       IF (IN24 .EQ. 24)  THEN
 129          IN24 = 0
 130          KOUNT = KOUNT + NSKIP
 131          DO 90 ISK= 1, NSKIP
 132          UNI = SEEDS(J24) - SEEDS(I24) - CARRY
 133          IF (UNI .LT. 0.)  THEN
 134             UNI = UNI + 1.0
 135             CARRY = TWOM24
 136          ELSE
 137             CARRY = 0.
 138          ENDIF
 139          SEEDS(I24) = UNI
 140          I24 = NEXT(I24)
 141          J24 = NEXT(J24)
 142    90    CONTINUE
 143       ENDIF
 144   100 CONTINUE
 145       KOUNT = KOUNT + LENV
 146       IF (KOUNT .GE. IGIGA)  THEN
 147          MKOUNT = MKOUNT + 1
 148          KOUNT = KOUNT - IGIGA
 149       ENDIF
 150       RETURN
 151 C
 152 C           Entry to input and float integer seeds from previous run
 153       ENTRY RLUXIN(ISDEXT)
 154          NOTYET = .FALSE.
 155          TWOM24 = 1.
 156          DO 195 I= 1, 24
 157          NEXT(I) = I-1
 158   195    TWOM24 = TWOM24 * 0.5
 159          NEXT(1) = 24
 160          TWOM12 = TWOM24 * 4096.
 161       WRITE(6,'(A)') ' FULL INITIALIZATION OF RANLUX WITH 25 INTEGERS:'
 162       WRITE(6,'(5X,5I12)') ISDEXT
 163       DO 200 I= 1, 24
 164       SEEDS(I) = REAL(ISDEXT(I))*TWOM24
 165   200 CONTINUE
 166       CARRY = 0.
 167       IF (ISDEXT(25) .LT. 0)  CARRY = TWOM24
 168       ISD = IABS(ISDEXT(25))
 169       I24 = MOD(ISD,100)
 170       ISD = ISD/100
 171       J24 = MOD(ISD,100)
 172       ISD = ISD/100
 173       IN24 = MOD(ISD,100)
 174       ISD = ISD/100
 175       LUXLEV = ISD
 176         IF (LUXLEV .LE. MAXLEV) THEN
 177           NSKIP = NDSKIP(LUXLEV)
 178           WRITE (6,'(A,I2)') ' RANLUX LUXURY LEVEL SET BY RLUXIN TO: ',
 179      +                         LUXLEV
 180         ELSE  IF (LUXLEV .GE. 24) THEN
 181           NSKIP = LUXLEV - 24
 182           WRITE (6,'(A,I5)') ' RANLUX P-VALUE SET BY RLUXIN TO:',LUXLEV
 183         ELSE
 184           NSKIP = NDSKIP(MAXLEV)
 185           WRITE (6,'(A,I5)') ' RANLUX ILLEGAL LUXURY RLUXIN: ',LUXLEV
 186           LUXLEV = MAXLEV
 187         ENDIF
 188       INSEED = -1
 189       RETURN
 190 C
 191 C                    Entry to ouput seeds as integers
 192       ENTRY RLUXUT(ISDEXT)
 193       DO 300 I= 1, 24
 194          ISDEXT(I) = INT(SEEDS(I)*TWOP12*TWOP12)
 195   300 CONTINUE
 196       ISDEXT(25) = I24 + 100*J24 + 10000*IN24 + 1000000*LUXLEV
 197       IF (CARRY .GT. 0.)  ISDEXT(25) = -ISDEXT(25)
 198       RETURN
 199 C
 200 C                    Entry to output the "convenient" restart point
 201       ENTRY RLUXAT(LOUT,INOUT,K1,K2)
 202       LOUT = LUXLEV
 203       INOUT = INSEED
 204       K1 = KOUNT
 205       K2 = MKOUNT
 206       RETURN
 207 C
 208 C                    Entry to initialize from one or three integers
 209       ENTRY RLUXGO(LUX,INS,K1,K2)
 210          IF (LUX .LT. 0) THEN
 211             LUXLEV = LXDFLT
 212          ELSE IF (LUX .LE. MAXLEV) THEN
 213             LUXLEV = LUX
 214          ELSE IF (LUX .LT. 24 .OR. LUX .GT. 2000) THEN
 215             LUXLEV = MAXLEV
 216             WRITE (6,'(A,I7)') ' RANLUX ILLEGAL LUXURY RLUXGO: ',LUX
 217          ELSE
 218             LUXLEV = LUX
 219             DO 310 ILX= 0, MAXLEV
 220               IF (LUX .EQ. NDSKIP(ILX)+24)  LUXLEV = ILX
 221   310       CONTINUE
 222          ENDIF
 223       IF (LUXLEV .LE. MAXLEV)  THEN
 224          NSKIP = NDSKIP(LUXLEV)
 225          WRITE(6,'(A,I2,A,I4)') ' RANLUX LUXURY LEVEL SET BY RLUXGO :',
 226      +        LUXLEV,'     P=', NSKIP+24
 227       ELSE
 228           NSKIP = LUXLEV - 24
 229           WRITE (6,'(A,I5)') ' RANLUX P-VALUE SET BY RLUXGO TO:',LUXLEV
 230       ENDIF
 231       IN24 = 0
 232       IF (INS .LT. 0)  WRITE (6,'(A)')
 233      +   ' Illegal initialization by RLUXGO, negative input seed'
 234       IF (INS .GT. 0)  THEN
 235         JSEED = INS
 236         WRITE(6,'(A,3I12)') ' RANLUX INITIALIZED BY RLUXGO FROM SEEDS',
 237      +      JSEED, K1,K2
 238       ELSE
 239         JSEED = JSDFLT
 240         WRITE(6,'(A)')' RANLUX INITIALIZED BY RLUXGO FROM DEFAULT SEED'
 241       ENDIF
 242       INSEED = JSEED
 243       NOTYET = .FALSE.
 244       TWOM24 = 1.
 245          DO 325 I= 1, 24
 246            TWOM24 = TWOM24 * 0.5
 247          K = JSEED/53668
 248          JSEED = 40014*(JSEED-K*53668) -K*12211
 249          IF (JSEED .LT. 0)  JSEED = JSEED+ICONS
 250          ISEEDS(I) = MOD(JSEED,ITWO24)
 251   325    CONTINUE
 252       TWOM12 = TWOM24 * 4096.
 253          DO 350 I= 1,24
 254          SEEDS(I) = REAL(ISEEDS(I))*TWOM24
 255          NEXT(I) = I-1
 256   350    CONTINUE
 257       NEXT(1) = 24
 258       I24 = 24
 259       J24 = 10
 260       CARRY = 0.
 261       IF (SEEDS(24) .EQ. 0.) CARRY = TWOM24
 262 C        If restarting at a break point, skip K1 + IGIGA*K2
 263 C        Note that this is the number of numbers delivered to
 264 C        the user PLUS the number skipped (if luxury .GT. 0).
 265       KOUNT = K1
 266       MKOUNT = K2
 267       IF (K1+K2 .NE. 0)  THEN
 268         DO 500 IOUTER= 1, K2+1
 269           INNER = IGIGA
 270           IF (IOUTER .EQ. K2+1)  INNER = K1
 271           DO 450 ISK= 1, INNER
 272             UNI = SEEDS(J24) - SEEDS(I24) - CARRY
 273             IF (UNI .LT. 0.)  THEN
 274                UNI = UNI + 1.0
 275                CARRY = TWOM24
 276             ELSE
 277                CARRY = 0.
 278             ENDIF
 279             SEEDS(I24) = UNI
 280             I24 = NEXT(I24)
 281             J24 = NEXT(J24)
 282   450     CONTINUE
 283   500   CONTINUE
 284 C         Get the right value of IN24 by direct calculation
 285         IN24 = MOD(KOUNT, NSKIP+24)
 286         IF (MKOUNT .GT. 0)  THEN
 287            IZIP = MOD(IGIGA, NSKIP+24)
 288            IZIP2 = MKOUNT*IZIP + IN24
 289            IN24 = MOD(IZIP2, NSKIP+24)
 290         ENDIF
 291 C       Now IN24 had better be between zero and 23 inclusive
 292         IF (IN24 .GT. 23) THEN
 293            WRITE (6,'(A/A,3I11,A,I5)')
 294      +    '  Error in RESTARTING with RLUXGO:','  The values', INS,
 295      +     K1, K2, ' cannot occur at luxury level', LUXLEV
 296            IN24 = 0
 297         ENDIF
 298       ENDIF
 299       RETURN
 300       END