MEVSIM/ranlux.F

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