tenstream solver
ranlux.F90
Go to the documentation of this file.
1 module m_ranlux
2 
3  !implicit none
4 
5  private
6  public :: ranlux, rluxgo
7 
8  dimension seeds(24), iseeds(24)
9  parameter(maxlev=4, lxdflt=3)
10  dimension ndskip(0:maxlev)
11  dimension next(24)
12  parameter(twop12=4096., igiga=1000000000, jsdflt=314159265)
13  parameter(itwo24=2**24, icons=2147483563)
14  save notyet, i24, j24, carry, seeds, twom24, twom12, luxlev
15  save nskip, ndskip, in24, next, kount, mkount, inseed
16  integer luxlev
17  logical notyet
18  data notyet, luxlev, in24, kount, mkount/.true., lxdflt, 0, 0, 0/
19  data i24, j24, carry/24, 10, 0./
20  ! default
21  ! Luxury Level 0 1 2 *3* 4
22  data ndskip/0, 24, 73, 199, 365/
23  !orresponds to p=24 48 97 223 389
24  ! time factor 1 2 3 6 10 on slow workstation
25  ! 1 1.5 2 3 5 on fast mainframe
26  !
27 contains
28  subroutine ranlux(RVEC, LENV)
29 ! Subtract-and-borrow random number generator proposed by
30 ! Marsaglia and Zaman, implemented by F. James with the name
31 ! RCARRY in 1991, and later improved by Martin Luescher
32 ! in 1993 to produce "Luxury Pseudorandom Numbers".
33 ! Fortran 77 coded by F. James, 1993
34 !
35 ! references:
36 ! M. Luscher, Computer Physics Communications 79 (1994) 100
37 ! F. James, Computer Physics Communications 79 (1994) 111
38 !
39 ! LUXURY LEVELS.
40 ! ------ ------ The available luxury levels are:
41 !
42 ! level 0 (p=24): equivalent to the original RCARRY of Marsaglia
43 ! and Zaman, very long period, but fails many tests.
44 ! level 1 (p=48): considerable improvement in quality over level 0,
45 ! now passes the gap test, but still fails spectral test.
46 ! level 2 (p=97): passes all known tests, but theoretically still
47 ! defective.
48 ! level 3 (p=223): DEFAULT VALUE. Any theoretically possible
49 ! correlations have very small chance of being observed.
50 ! level 4 (p=389): highest possible luxury, all 24 bits chaotic.
51 !
52 !!!! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
53 !!!! Calling sequences for RANLUX: ++
54 !!!! CALL RANLUX (RVEC, LEN) returns a vector RVEC of LEN ++
55 !!!! 32-bit random floating point numbers between ++
56 !!!! zero (not included) and one (also not incl.). ++
57 !!!! CALL RLUXGO(LUX,INT,K1,K2) initializes the generator from ++
58 !!!! one 32-bit integer INT and sets Luxury Level LUX ++
59 !!!! which is integer between zero and MAXLEV, or if ++
60 !!!! LUX .GT. 24, it sets p=LUX directly. K1 and K2 ++
61 !!!! should be set to zero unless restarting at a break++
62 !!!! point given by output of RLUXAT (see RLUXAT). ++
63 !!!! CALL RLUXAT(LUX,INT,K1,K2) gets the values of four integers++
64 !!!! which can be used to restart the RANLUX generator ++
65 !!!! at the current point by calling RLUXGO. K1 and K2++
66 !!!! specify how many numbers were generated since the ++
67 !!!! initialization with LUX and INT. The restarting ++
68 !!!! skips over K1+K2*E9 numbers, so it can be long.++
69 !!!! A more efficient but less convenient way of restarting is by: ++
70 !!!! CALL RLUXIN(ISVEC) restarts the generator from vector ++
71 !!!! ISVEC of 25 32-bit integers (see RLUXUT) ++
72 !!!! CALL RLUXUT(ISVEC) outputs the current values of the 25 ++
73 !!!! 32-bit integer seeds, to be used for restarting ++
74 !!!! ISVEC must be dimensioned 25 in the calling program ++
75 !!!! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
76  dimension rvec(lenv)
77 ! NOTYET is .TRUE. if no initialization has been performed yet.
78 ! Default Initialization by Multiplicative Congruential
79  if (notyet) then
80  notyet = .false.
81  jseed = jsdflt
82  inseed = jseed
83  write (6, '(A,I12)') ' RANLUX DEFAULT INITIALIZATION: ', jseed
84  luxlev = lxdflt
85  nskip = ndskip(luxlev)
86  lp = nskip + 24
87  in24 = 0
88  kount = 0
89  mkount = 0
90  write (6, '(A,I2,A,I4)') ' RANLUX DEFAULT LUXURY LEVEL = ', &
91  +luxlev, ' p =', lp
92  twom24 = 1.
93  do 25 i = 1, 24
94  twom24 = twom24 * 0.5
95  k = jseed / 53668
96  jseed = 40014 * (jseed - k * 53668) - k * 12211
97  if (jseed .lt. 0) jseed = jseed + icons
98  iseeds(i) = mod(jseed, itwo24)
99 25 continue
100  twom12 = twom24 * 4096.
101  do 50 i = 1, 24
102  seeds(i) = real(ISEEDS(I)) * TWOM24
103  next(i) = i - 1
104 50 continue
105  next(1) = 24
106  i24 = 24
107  j24 = 10
108  carry = 0.
109  if (abs(seeds(24)) .lt. tiny(seeds)) carry = twom24
110  end if
111 !
112 ! The Generator proper: "Subtract-with-borrow",
113 ! as proposed by Marsaglia and Zaman,
114 ! Florida State University, March, 1989
115 !
116  do 100 ivec = 1, lenv
117  uni = seeds(j24) - seeds(i24) - carry
118  if (uni .lt. 0.) then
119  uni = uni + 1.0
120  carry = twom24
121  else
122  carry = 0.
123  end if
124  seeds(i24) = uni
125  i24 = next(i24)
126  j24 = next(j24)
127  rvec(ivec) = uni
128 ! small numbers (with less than 12 "significant" bits) are "padded".
129  if (uni .lt. twom12) then
130  rvec(ivec) = rvec(ivec) + twom24 * seeds(j24)
131 ! and zero is forbidden in case someone takes a logarithm
132  if (abs(rvec(ivec)) .lt. tiny(rvec)) rvec(ivec) = twom24 * twom24
133  end if
134 ! Skipping to luxury. As proposed by Martin Luscher.
135  in24 = in24 + 1
136  if (in24 .eq. 24) then
137  in24 = 0
138  kount = kount + nskip
139  do 90 isk = 1, nskip
140  uni = seeds(j24) - seeds(i24) - carry
141  if (uni .lt. 0.) then
142  uni = uni + 1.0
143  carry = twom24
144  else
145  carry = 0.
146  end if
147  seeds(i24) = uni
148  i24 = next(i24)
149  j24 = next(j24)
150 90 continue
151  end if
152 100 continue
153  kount = kount + lenv
154  if (kount .ge. igiga) then
155  mkount = mkount + 1
156  kount = kount - igiga
157  end if
158  return
159  end subroutine ranlux
160 !
161 ! Entry to input and float integer seeds from previous run
162  subroutine rluxin(ISDEXT)
163  dimension isdext(25)
164 !--------------------------added to fix `bug'---------------------------
165  if (notyet) then
166 !*****$#--1----_----2----_----3----_----4----_----5----_----6----_----7--
167 ! WRITE(6,'(A)') ' PROPER RESULTS ONLY WITH INITIALISATION FROM
168 ! $25 INTEGERS OBTAINED WITH RLUXUT'
169  notyet = .false.
170 !-----------------------------------------------------------------------
171  end if
172  twom24 = 1.
173  do i = 1, 24
174  next(i) = i - 1
175  twom24 = twom24 * 0.5
176  end do
177  next(1) = 24
178  twom12 = twom24 * 4096.
179  write (6, '(A)') ' FULL INITIALIZATION OF RANLUX WITH 25 INTEGERS:'
180  write (6, '(5X,5I12)') isdext
181  do 200 i = 1, 24
182  seeds(i) = real(ISDEXT(I)) * TWOM24
183 200 continue
184  carry = 0.
185  if (isdext(25) .lt. 0) carry = twom24
186  isd = iabs(isdext(25))
187  i24 = mod(isd, 100)
188  isd = isd / 100
189  j24 = mod(isd, 100)
190  isd = isd / 100
191  in24 = mod(isd, 100)
192  isd = isd / 100
193  luxlev = isd
194  if (luxlev .le. maxlev) then
195  nskip = ndskip(luxlev)
196  write (6, '(A,I2)') ' RANLUX LUXURY LEVEL SET BY RLUXIN TO: ', &
197  +luxlev
198  else if (luxlev .ge. 24) then
199  nskip = luxlev - 24
200  write (6, '(A,I5)') ' RANLUX P-VALUE SET BY RLUXIN TO:', luxlev
201  else
202  nskip = ndskip(maxlev)
203  write (6, '(A,I5)') ' RANLUX ILLEGAL LUXURY RLUXIN: ', luxlev
204  luxlev = maxlev
205  end if
206  inseed = -1
207  end subroutine rluxin
208 !
209 ! Entry to ouput seeds as integers
210  subroutine rluxut(ISDEXT)
211  dimension isdext(25)
212  do i = 1, 24
213  isdext(i) = int(seeds(i) * twop12 * twop12)
214  end do
215  isdext(25) = i24 + 100 * j24 + 10000 * in24 + 1000000 * luxlev
216  if (carry .gt. 0.) isdext(25) = -isdext(25)
217  end subroutine rluxut
218 !
219 ! Entry to output the "convenient" restart point
220  subroutine rluxat(LOUT, INOUT, K1, K2)
221  lout = luxlev
222  inout = inseed
223  k1 = kount
224  k2 = mkount
225  end subroutine rluxat
226 !
227 ! Entry to initialize from one or three integers
228  subroutine rluxgo(LUX, INS, K1, K2)
229  if (lux .lt. 0) then
230  luxlev = lxdflt
231  else if (lux .le. maxlev) then
232  luxlev = lux
233  else if (lux .lt. 24 .or. lux .gt. 2000) then
234  luxlev = maxlev
235  write (6, '(A,I7)') ' RANLUX ILLEGAL LUXURY RLUXGO: ', lux
236  else
237  luxlev = lux
238  do 310 ilx = 0, maxlev
239  if (lux .eq. ndskip(ilx) + 24) luxlev = ilx
240 310 continue
241  end if
242  if (luxlev .le. maxlev) then
243  nskip = ndskip(luxlev)
244 ! WRITE(6,'(A,I2,A,I4)') ' RANLUX LUXURY LEVEL SET BY RLUXGO :',
245 ! + LUXLEV,' P=', NSKIP+24
246  else
247  nskip = luxlev - 24
248  write (6, '(A,I5)') ' RANLUX P-VALUE SET BY RLUXGO TO:', luxlev
249  end if
250  in24 = 0
251  if (ins .lt. 0) write (6, '(A)') ' Illegal initialization by RLUXGO, negative input seed'
252  if (ins .gt. 0) then
253  jseed = ins
254 ! WRITE(6,'(A,3I12)') ' RANLUX INITIALIZED BY RLUXGO FROM SEEDS',
255 ! + JSEED, K1,K2
256  else
257  jseed = jsdflt
258  write (6, '(A)') ' RANLUX INITIALIZED BY RLUXGO FROM DEFAULT SEED'
259  end if
260  inseed = jseed
261  notyet = .false.
262  twom24 = 1.
263  do 325 i = 1, 24
264  twom24 = twom24 * 0.5
265  k = jseed / 53668
266  jseed = 40014 * (jseed - k * 53668) - k * 12211
267  if (jseed .lt. 0) jseed = jseed + icons
268  iseeds(i) = mod(jseed, itwo24)
269 325 continue
270  twom12 = twom24 * 4096.
271  do 350 i = 1, 24
272  seeds(i) = real(ISEEDS(I)) * TWOM24
273  next(i) = i - 1
274 350 continue
275  next(1) = 24
276  i24 = 24
277  j24 = 10
278  carry = 0.
279  if (abs(seeds(24)) .lt. tiny(seeds)) carry = twom24
280 ! If restarting at a break point, skip K1 + IGIGA*K2
281 ! Note that this is the number of numbers delivered to
282 ! the user PLUS the number skipped (if luxury .GT. 0).
283  kount = k1
284  mkount = k2
285  if (k1 + k2 .ne. 0) then
286  do 500 iouter = 1, k2 + 1
287  inner = igiga
288  if (iouter .eq. k2 + 1) inner = k1
289  do 450 isk = 1, inner
290  uni = seeds(j24) - seeds(i24) - carry
291  if (uni .lt. 0.) then
292  uni = uni + 1.0
293  carry = twom24
294  else
295  carry = 0.
296  end if
297  seeds(i24) = uni
298  i24 = next(i24)
299  j24 = next(j24)
300 450 continue
301 500 continue
302 ! Get the right value of IN24 by direct calculation
303  in24 = mod(kount, nskip + 24)
304  if (mkount .gt. 0) then
305  izip = mod(igiga, nskip + 24)
306  izip2 = mkount * izip + in24
307  in24 = mod(izip2, nskip + 24)
308  end if
309 ! Now IN24 had better be between zero and 23 inclusive
310  if (in24 .gt. 23) then
311  write (6, '(A/A,3I11,A,I5)') ' Error in RESTARTING with RLUXGO:',&
312  & ' The values', ins, k1, k2, ' cannot occur at luxury level', luxlev
313  in24 = 0
314  end if
315  end if
316  end subroutine rluxgo
317  end module
m_ranlux::notyet
logical, save notyet
Definition: ranlux.F90:17
m_ranlux::rluxgo
subroutine, public rluxgo(LUX, INS, K1, K2)
Definition: ranlux.F90:229
m_ranlux::rluxat
subroutine rluxat(LOUT, INOUT, K1, K2)
Definition: ranlux.F90:221
m_ranlux::rluxut
subroutine rluxut(ISDEXT)
Definition: ranlux.F90:211
m_ranlux
Definition: ranlux.F90:1
m_ranlux::luxlev
integer, save luxlev
Definition: ranlux.F90:16
m_ranlux::ranlux
subroutine, public ranlux(RVEC, LENV)
Definition: ranlux.F90:29
m_ranlux::rluxin
subroutine rluxin(ISDEXT)
Definition: ranlux.F90:163