source: branches/jerome/src_flexwrf_v3.1/ranlux.f90 @ 16

Last change on this file since 16 was 16, checked in by jebri, 11 years ago

sources for flexwrf v3.1
File size: 16.9 KB
Line 
1 
2    MODULE luxury
3
4!     Subtract-and-borrow random number generator proposed by
5!     Marsaglia and Zaman, implemented by F. James with the name
6!     RCARRY in 1991, and later improved by Martin Luescher
7!     in 1993 to produce "Luxury Pseudorandom Numbers".
8!     Fortran 77 coded by F. James, 1993
9
10!  References:
11!  M. Luscher, Computer Physics Communications  79 (1994) 100
12!  F. James, Computer Physics Communications 79 (1994) 111
13
14!   LUXURY LEVELS.
15!   ------ ------      The available luxury levels are:
16
17!  level 0  (p=24): equivalent to the original RCARRY of Marsaglia
18!           and Zaman, very long period, but fails many tests.
19!  level 1  (p=48): considerable improvement in quality over level 0,
20!           now passes the gap test, but still fails spectral test.
21!  level 2  (p=97): passes all known tests, but theoretically still
22!           defective.
23!  level 3  (p=223): DEFAULT VALUE.  Any theoretically possible
24!           correlations have very small chance of being observed.
25!  level 4  (p=389): highest possible luxury, all 24 bits chaotic.
26
27!!!! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
28!!!!  Calling sequences for RANLUX:                                  ++
29!!!!      CALL RANLUX (RVEC, LEN)   returns a vector RVEC of LEN     ++
30!!!!                   32-bit random floating point numbers between  ++
31!!!!                   zero (not included) and one (also not incl.). ++
32!!!!      CALL RLUXGO(LUX,INT,K1,K2) initializes the generator from  ++
33!!!!               one 32-bit integer INT and sets Luxury Level LUX  ++
34!!!!               which is integer between zero and MAXLEV, or if   ++
35!!!!               LUX .GT. 24, it sets p=LUX directly.  K1 and K2   ++
36!!!!               should be set to zero unless restarting at a break++
37!!!!               point given by output of RLUXAT (see RLUXAT).     ++
38!!!!      CALL RLUXAT(LUX,INT,K1,K2) gets the values of four integers++
39!!!!               which can be used to restart the RANLUX generator ++
40!!!!               at the current point by calling RLUXGO.  K1 and K2++
41!!!!               specify how many numbers were generated since the ++
42!!!!               initialization with LUX and INT.  The restarting  ++
43!!!!               skips over  K1+K2*E9   numbers, so it can be long.++
44!!!!   A more efficient but less convenient way of restarting is by: ++
45!!!!      CALL RLUXIN(ISVEC)    restarts the generator from vector   ++
46!!!!                   ISVEC of 25 32-bit integers (see RLUXUT)      ++
47!!!!      CALL RLUXUT(ISVEC)    outputs the current values of the 25 ++
48!!!!                 32-bit integer seeds, to be used for restarting ++
49!!!!      ISVEC must be dimensioned 25 in the calling program        ++
50!!!! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
51
52IMPLICIT NONE
53
54INTEGER            :: iseeds(24), isdext(25)
55INTEGER, PARAMETER :: maxlev = 4, lxdflt = 3
56INTEGER            :: ndskip(0:maxlev) = (/ 0, 24, 73, 199, 365 /)
57INTEGER            :: next(24), igiga = 1000000000, jsdflt = 314159265
58REAL, PARAMETER    :: twop12 = 4096.
59INTEGER, PARAMETER :: itwo24 = 2**24, icons = 2147483563
60INTEGER            :: luxlev = lxdflt, nskip, inseed, jseed
61LOGICAL            :: notyet = .true.
62INTEGER            :: in24 = 0, kount = 0, mkount = 0, i24 = 24, j24 = 10
63REAL               :: seeds(24), carry = 0., twom24, twom12
64
65!                            default
66!  Luxury Level     0   1   2  *3*    4
67!    ndskip        /0, 24, 73, 199, 365/
68! Corresponds to p=24  48  97  223  389
69!     time factor   1   2   3    6   10   on slow workstation
70!                   1 1.5   2    3    5   on fast mainframe
71
72PUBLIC notyet, i24, j24, carry, seeds, twom24, twom12, luxlev
73PUBLIC nskip, ndskip, in24, next, kount, mkount, inseed
74
75
76CONTAINS
77
78
79SUBROUTINE ranlux(rvec, lenv)
80
81IMPLICIT NONE
82
83INTEGER, INTENT(IN) :: lenv
84REAL, INTENT(OUT)   :: rvec(lenv)
85
86!     Local variables
87
88INTEGER             :: i, k, lp, ivec, isk
89REAL                :: uni
90
91!  NOTYET is .TRUE. if no initialization has been performed yet.
92!              Default Initialization by Multiplicative Congruential
93
94IF (notyet) THEN
95  notyet = .false.
96  jseed = jsdflt
97  inseed = jseed
98  WRITE (6,'(A,I12)') ' RANLUX DEFAULT INITIALIZATION: ', jseed
99  luxlev = lxdflt
100  nskip = ndskip(luxlev)
101  lp = nskip + 24
102  in24 = 0
103  kount = 0
104  mkount = 0
105  WRITE (6,'(A,I2,A,I4)') ' RANLUX DEFAULT LUXURY LEVEL =  ', luxlev,   &
106                          '    p =', lp
107  twom24 = 1.
108  DO i = 1, 24
109    twom24 = twom24 * 0.5
110    k = jseed / 53668
111    jseed = 40014 * (jseed-k*53668) - k * 12211
112    IF (jseed.LT.0) jseed = jseed + icons
113    iseeds(i) = MOD(jseed,itwo24)
114  END DO
115  twom12 = twom24 * 4096.
116  DO i = 1, 24
117    seeds(i) = REAL(iseeds(i)) * twom24
118    next(i) = i - 1
119  END DO
120  next(1) = 24
121  i24 = 24
122  j24 = 10
123  carry = 0.
124  IF (seeds(24).EQ.0.) carry = twom24
125END IF
126
127!          The Generator proper: "Subtract-with-borrow",
128!          as proposed by Marsaglia and Zaman,
129!          Florida State University, March, 1989
130
131DO ivec = 1, lenv
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  END IF
139  seeds(i24) = uni
140  i24 = next(i24)
141  j24 = next(j24)
142  rvec(ivec) = uni
143!  small numbers (with less than 12 "significant" bits) are "padded".
144  IF (uni.LT.twom12) THEN
145    rvec(ivec) = rvec(ivec) + twom24 * seeds(j24)
146!        and zero is forbidden in case someone takes a logarithm
147    IF (rvec(ivec).EQ.0.) rvec(ivec) = twom24 * twom24
148  END IF
149!        Skipping to luxury.  As proposed by Martin Luscher.
150  in24 = in24 + 1
151  IF (in24.EQ.24) THEN
152    in24 = 0
153    kount = kount + nskip
154    DO isk = 1, nskip
155      uni = seeds(j24) - seeds(i24) - carry
156      IF (uni.LT.0.) THEN
157        uni = uni + 1.0
158        carry = twom24
159      ELSE
160        carry = 0.
161      END IF
162      seeds(i24) = uni
163      i24 = next(i24)
164      j24 = next(j24)
165    END DO
166  END IF
167END DO
168kount = kount + lenv
169IF (kount.GE.igiga) THEN
170  mkount = mkount + 1
171  kount = kount - igiga
172END IF
173RETURN
174
175END SUBROUTINE ranlux
176
177
178!           Subroutine to input and float integer seeds from previous run
179SUBROUTINE rluxin
180!     the following IF BLOCK added by Phillip Helbig, based on conversation
181!     with Fred James; an equivalent correction has been published by James.
182
183IMPLICIT NONE
184
185!     Local variables
186
187INTEGER             :: i, isd
188
189IF (notyet) THEN
190  WRITE (6,'(A)') ' Proper results ONLY with initialisation from 25 ',  &
191  'integers obtained with RLUXUT'
192  notyet = .false.
193END IF
194
195twom24 = 1.
196DO i = 1, 24
197  next(i) = i - 1
198  twom24 = twom24 * 0.5
199END DO
200next(1) = 24
201twom12 = twom24 * 4096.
202WRITE (6,'(A)') ' FULL INITIALIZATION OF RANLUX WITH 25 INTEGERS:'
203WRITE (6,'(5X,5I12)') isdext
204DO i = 1, 24
205  seeds(i) = REAL(isdext(i)) * twom24
206END DO
207carry = 0.
208IF (isdext(25).LT.0) carry = twom24
209isd = IABS(isdext(25))
210i24 = MOD(isd,100)
211isd = isd / 100
212j24 = MOD(isd,100)
213isd = isd / 100
214in24 = MOD(isd,100)
215isd = isd / 100
216luxlev = isd
217IF (luxlev.LE.maxlev) THEN
218  nskip = ndskip(luxlev)
219  WRITE (6,'(A,I2)') ' RANLUX LUXURY LEVEL SET BY RLUXIN TO: ', luxlev
220ELSE IF (luxlev.GE.24) THEN
221  nskip = luxlev - 24
222  WRITE (6,'(A,I5)') ' RANLUX P-VALUE SET BY RLUXIN TO:', luxlev
223ELSE
224  nskip = ndskip(maxlev)
225  WRITE (6,'(A,I5)') ' RANLUX ILLEGAL LUXURY RLUXIN: ', luxlev
226  luxlev = maxlev
227END IF
228inseed = -1
229RETURN
230
231END SUBROUTINE rluxin
232
233
234!                    Subroutine to ouput seeds as integers
235SUBROUTINE rluxut
236
237IMPLICIT NONE
238
239!     Local variables
240
241INTEGER             :: i
242
243DO i = 1, 24
244  isdext(i) = INT(seeds(i)*twop12*twop12)
245END DO
246isdext(25) = i24 + 100 * j24 + 10000 * in24 + 1000000 * luxlev
247IF (carry.GT.0.) isdext(25) = -isdext(25)
248RETURN
249
250END SUBROUTINE rluxut
251
252
253!                    Subroutine to output the "convenient" restart point
254SUBROUTINE rluxat(lout, inout, k1, k2)
255
256IMPLICIT NONE
257
258INTEGER, INTENT(OUT) :: lout, inout, k1, k2
259
260lout = luxlev
261inout = inseed
262k1 = kount
263k2 = mkount
264RETURN
265
266END SUBROUTINE rluxat
267
268
269!                    Subroutine to initialize from one or three integers
270SUBROUTINE rluxgo(lux, ins, k1, k2)
271
272IMPLICIT NONE
273
274INTEGER, INTENT(IN) :: lux, ins, k1, k2
275
276!     Local variables
277
278INTEGER             :: ilx, i, iouter, isk, k, inner, izip, izip2
279REAL                :: uni
280
281IF (lux.LT.0) THEN
282  luxlev = lxdflt
283ELSE IF (lux.LE.maxlev) THEN
284  luxlev = lux
285ELSE IF (lux.LT.24.OR.lux.GT.2000) THEN
286  luxlev = maxlev
287  WRITE (6,'(A,I7)') ' RANLUX ILLEGAL LUXURY RLUXGO: ', lux
288ELSE
289  luxlev = lux
290  DO ilx = 0, maxlev
291    IF (lux.EQ.ndskip(ilx)+24) luxlev = ilx
292  END DO
293END IF
294IF (luxlev.LE.maxlev) THEN
295  nskip = ndskip(luxlev)
296  WRITE (6,'(A,I2,A,I4)') ' RANLUX LUXURY LEVEL SET BY RLUXGO :', luxlev,  &
297                          '     P=', nskip + 24
298ELSE
299  nskip = luxlev - 24
300  WRITE (6,'(A,I5)') ' RANLUX P-VALUE SET BY RLUXGO TO:', luxlev
301END IF
302in24 = 0
303IF (ins.LT.0) WRITE (6,'(A)') &
304              ' Illegal initialization by RLUXGO, negative input seed'
305IF (ins.GT.0) THEN
306  jseed = ins
307  WRITE (6,'(A,3I12)') ' RANLUX INITIALIZED BY RLUXGO FROM SEEDS', jseed, k1, k2
308ELSE
309  jseed = jsdflt
310  WRITE (6,'(A)') ' RANLUX INITIALIZED BY RLUXGO FROM DEFAULT SEED'
311END IF
312inseed = jseed
313notyet = .false.
314twom24 = 1.
315DO i = 1, 24
316  twom24 = twom24 * 0.5
317  k = jseed / 53668
318  jseed = 40014 * (jseed-k*53668) - k * 12211
319  IF (jseed.LT.0) jseed = jseed + icons
320  iseeds(i) = MOD(jseed,itwo24)
321END DO
322twom12 = twom24 * 4096.
323DO i = 1, 24
324  seeds(i) = REAL(iseeds(i)) * twom24
325  next(i) = i - 1
326END DO
327next(1) = 24
328i24 = 24
329j24 = 10
330carry = 0.
331IF (seeds(24).EQ.0.) carry = twom24
332!        If restarting at a break point, skip K1 + IGIGA*K2
333!        Note that this is the number of numbers delivered to
334!        the user PLUS the number skipped (if luxury .GT. 0).
335kount = k1
336mkount = k2
337IF (k1+k2.NE.0) THEN
338  DO iouter = 1, k2 + 1
339    inner = igiga
340    IF (iouter.EQ.k2+1) inner = k1
341    DO isk = 1, inner
342      uni = seeds(j24) - seeds(i24) - carry
343      IF (uni.LT.0.) THEN
344        uni = uni + 1.0
345        carry = twom24
346      ELSE
347        carry = 0.
348      END IF
349      seeds(i24) = uni
350      i24 = next(i24)
351      j24 = next(j24)
352    END DO
353  END DO
354!         Get the right value of IN24 by direct calculation
355  in24 = MOD(kount,nskip+24)
356  IF (mkount.GT.0) THEN
357    izip = MOD(igiga, nskip+24)
358    izip2 = mkount * izip + in24
359    in24 = MOD(izip2, nskip+24)
360  END IF
361!       Now IN24 had better be between zero and 23 inclusive
362  IF (in24.GT.23) THEN
363    WRITE (6,'(A/A,3I11,A,I5)') &
364               '  Error in RESTARTING with RLUXGO:', '  The values', ins, &
365               k1, k2, ' cannot occur at luxury level', luxlev
366    in24 = 0
367  END IF
368END IF
369RETURN
370
371END SUBROUTINE rluxgo
372
373
374END MODULE luxury
375
376
377
378!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
379
380subroutine luxtst
381!         Exercise for the RANLUX Pseudorandom number generator.
382
383USE luxury
384
385IMPLICIT NONE
386
387REAL    :: rvec(1000)
388INTEGER :: i1, i2, i3, i4, li
389
390!         check that we get the right numbers (machine-indep.)
391WRITE (6,'(/A)') '  CALL RANLUX(RVEC,100)'
392CALL ranlux(rvec,100)
393WRITE (6,'(A/9X,5F12.8)') ' RANLUX default numbers   1-  5:', rvec(1:5)
394CALL ranlux(rvec,100)
395WRITE (6,'(A/9X,5F12.8)') ' RANLUX default numbers 101-105:', rvec(1:5)
396
397WRITE (6,'(/A)') ' CALL RLUXGO(0,0,0,0)'
398CALL rluxgo(0,0,0,0)
399CALL ranlux(rvec,100)
400WRITE (6,'(A/9X,5F12.8)') ' RANLUX luxury level 0,   1-  5:', rvec(1:5)
401CALL ranlux(rvec,100)
402WRITE (6,'(A/9X,5F12.8)') ' RANLUX luxury level 0, 101-105:', rvec(1:5)
403
404WRITE (6,'(/A)') '   CALL RLUXGO(389,1,0,0)'
405CALL rluxgo(389,1,0,0)
406CALL ranlux(rvec,100)
407WRITE (6,'(A/9X,5F12.8)') ' RANLUX luxury p=389,   1-  5:', rvec(1:5)
408CALL ranlux(rvec,100)
409WRITE (6,'(A/9X,5F12.8)') ' RANLUX luxury p=389, 101-105:', rvec(1:5)
410
411WRITE (6,'(/A)') '  CALL RLUXGO(75,0,0,0)'
412CALL rluxgo(75,0,0,0)
413CALL ranlux(rvec,100)
414WRITE (6,'(A/9X,5F12.8)') ' RANLUX luxury p= 75,   1-  5:', rvec(1:5)
415CALL ranlux(rvec,100)
416WRITE (6,'(A/9X,5F12.8)') ' RANLUX luxury p= 75, 101-105:', rvec(1:5)
417
418WRITE (6,'(/A)') '  test restarting from the full vector'
419CALL rluxut
420WRITE (6,'(/A/(1X,5I14))') '  current RANLUX status saved:', isdext
421CALL ranlux(rvec,100)
422WRITE (6,'(A/9X,5F12.8)') ' RANLUX numbers 1- 5:', rvec(1:5)
423CALL ranlux(rvec,100)
424WRITE (6,'(A/9X,5F12.8)') ' RANLUX numbers 101-105:', rvec(1:5)
425
426WRITE (6,'(/A)') '   previous RANLUX status will be restored'
427CALL rluxin
428CALL ranlux(rvec,100)
429WRITE (6,'(A/9X,5F12.8)') ' RANLUX numbers 1- 5:', rvec(1:5)
430CALL ranlux(rvec,100)
431WRITE (6,'(A/9X,5F12.8)') ' RANLUX numbers 101-105:', rvec(1:5)
432
433WRITE (6,'(/A)') '     test the restarting by skipping'
434CALL rluxgo(4,7674985,0,0)
435CALL rluxat(i1,i2,i3,i4)
436WRITE (6,'(A,4I10)') '  RLUXAT values =', i1, i2, i3, i4
437DO li = 1, 10
438  CALL ranlux(rvec,1000)
439END DO
440CALL rluxat(i1,i2,i3,i4)
441WRITE (6,'(A,4I10)') '  RLUXAT values =', i1, i2, i3, i4
442CALL ranlux(rvec,200)
443WRITE (6,'(A,2F10.6)') '  Next and 200th numbers are:', rvec(1), rvec(200)
444CALL rluxgo(i1,i2,i3,i4)
445CALL ranlux(rvec,200)
446WRITE (6,'(A,2F10.6)') '  Next and 200th numbers are:', rvec(1), rvec(200)
447
448WRITE (6,'(/A)') ' The following should provoke an error message'
449CALL rluxgo(4,11111,31,0)
450STOP
451
452!   OUTPUT FROM THE ABOVE TEST PROGRAM SHOULD BE:
453!   --------------------------------------------
454!  CALL RANLUX(RVEC,100)
455! RANLUX DEFAULT INITIALIZATION:    314159265
456! RANLUX DEFAULT LUXURY LEVEL =   3      p = 223
457! RANLUX default numbers   1-  5:
458!           0.53981817  0.76155043  0.06029940  0.79600263  0.30631220
459! RANLUX default numbers 101-105:
460!           0.43156743  0.03774416  0.24897110  0.00147784  0.90274453
461
462!  CALL RLUXGO(0,0,0,0)
463! RANLUX LUXURY LEVEL SET BY RLUXGO : 0     P=  24
464! RANLUX INITIALIZED BY RLUXGO FROM DEFAULT SEED
465! RANLUX luxury level 0,   1-  5:
466!           0.53981817  0.76155043  0.06029940  0.79600263  0.30631220
467! RANLUX luxury level 0, 101-105:
468!           0.41538775  0.05330932  0.58195311  0.91397446  0.67034441
469
470!   CALL RLUXGO(389,1,0,0)
471! RANLUX LUXURY LEVEL SET BY RLUXGO : 4     P= 389
472! RANLUX INITIALIZED BY RLUXGO FROM SEEDS           1           0           0
473! RANLUX luxury p=389,   1-  5:
474!           0.94589490  0.47347850  0.95152789  0.42971975  0.09127384
475! RANLUX luxury p=389, 101-105:
476!           0.02618265  0.03775346  0.97274780  0.13302165  0.43126065
477
478!  CALL RLUXGO(75,0,0,0)
479! RANLUX P-VALUE SET BY RLUXGO TO:   75
480! RANLUX INITIALIZED BY RLUXGO FROM DEFAULT SEED
481! RANLUX luxury p= 75,   1-  5:
482!           0.53981817  0.76155043  0.06029940  0.79600263  0.30631220
483! RANLUX luxury p= 75, 101-105:
484!           0.25600731  0.23443210  0.59164381  0.59035838  0.07011414
485
486!  test restarting from the full vector
487
488!  current RANLUX status saved:
489!       16156027      16534309      15243811       2751687       6002207
490!        7979506       1301976       4567313       4305996       5872599
491!       12003090       2146823      12606367       4111505       5979640
492!       12739666      10489318      14036909      11729352       8061448
493!        7832659       6069758       3197719       1832730      75080216
494! RANLUX numbers 1- 5:
495!           0.22617835  0.60655993  0.86417443  0.43920082  0.23382509
496! RANLUX numbers 101-105:
497!           0.08107197  0.21466845  0.84856731  0.94078046  0.85626233
498
499!   previous RANLUX status will be restored
500! FULL INITIALIZATION OF RANLUX WITH 25 INTEGERS:
501!         16156027    16534309    15243811     2751687     6002207
502!          7979506     1301976     4567313     4305996     5872599
503!         12003090     2146823    12606367     4111505     5979640
504!         12739666    10489318    14036909    11729352     8061448
505!          7832659     6069758     3197719     1832730    75080216
506! RANLUX P-VALUE SET BY RLUXIN TO:   75
507! RANLUX numbers 1- 5:
508!           0.22617835  0.60655993  0.86417443  0.43920082  0.23382509
509! RANLUX numbers 101-105:
510!           0.08107197  0.21466845  0.84856731  0.94078046  0.85626233
511
512!     test the restarting by skipping
513! RANLUX LUXURY LEVEL SET BY RLUXGO : 4     P= 389
514! RANLUX INITIALIZED BY RLUXGO FROM SEEDS     7674985           0           0
515!  RLUXAT values =         4   7674985         0         0
516!  RLUXAT values =         4   7674985    161840         0
517!  Next and 200th numbers are:  0.019648  0.590586
518! RANLUX LUXURY LEVEL SET BY RLUXGO : 4     P= 389
519! RANLUX INITIALIZED BY RLUXGO FROM SEEDS     7674985      161840           0
520!  Next and 200th numbers are:  0.019648  0.590586
521
522! The following should provoke an error message
523! RANLUX LUXURY LEVEL SET BY RLUXGO : 4     P= 389
524! RANLUX INITIALIZED BY RLUXGO FROM SEEDS       11111          31           0
525!  Error in RESTARTING with RLUXGO:
526!  The values      11111         31          0 cannot occur at luxury level    4
527END subroutine luxtst
528
Note: See TracBrowser for help on using the repository browser.
hosted by ZAMG