Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Parallel (pseudo)random number generator (RNG) for multiple streams
10 : !> and substreams of random numbers.
11 : !>
12 : !> In detail, this RNG provides 2**64 random number streams each with a
13 : !> length of 2**127 resulting in a length of 2**191 for the total RNG.
14 : !> Moreover, each stream is divided in 2**51 substream of length 2**76.
15 : !> The stream lengths refer to the default precision of 32 bit random
16 : !> number, but also an extended precision of 53 bit per random number
17 : !> can be requested. In this case, two 32 bit random numbers are used
18 : !> to generate a 53 bit random number and therefore the stream length
19 : !> is halved when extended precision are requested.
20 : !>
21 : !> Usage hint:
22 : !>
23 : !> type(rng_stream_type) :: rng_stream
24 : !> rng_stream = rng_stream_type(name, ..., error=error)
25 : !>
26 : !> to generate the first stream. Optionally, you may define a different
27 : !> seed or create a stream of extended precision (53 bits). Then
28 : !>
29 : !> type(rng_stream_type) :: next_rng_stream
30 : !> next_rng_stream = rng_stream_type(name, last_rng_stream=rng_stream)
31 : !>
32 : !> to create all the following RNG streams w.r.t. the previous stream.
33 : !> The command line
34 : !>
35 : !> x = rng_stream%next(error=error)
36 : !>
37 : !> will provide the next real random number x between 0 and 1 and
38 : !>
39 : !> ix = rng_stream%next(low, high, error=error)
40 : !>
41 : !> the next integer random number ix between low and high from stream
42 : !> rng_stream. The default distribution type is a uniform distribution
43 : !> [0,1], but also other distribution types are available (see below).
44 : !>
45 : !> \par Literature
46 : !> P. L'Ecuyer, R. Simard, E. J. Chen, and W. D. Kelton,
47 : !> "An object-oriented random-number package with many long streams
48 : !> and substreams", Operations Research 50(6), 1073-1075 (2002)
49 : !> \author C++ code converted to Fortran 90/95 (18.05.2005, Matthias Krack)
50 : ! **************************************************************************************************
51 : MODULE parallel_rng_types
52 :
53 : USE kinds, ONLY: default_string_length,&
54 : dp
55 : USE string_utilities, ONLY: compress
56 : #include "../base/base_uses.f90"
57 :
58 : IMPLICIT NONE
59 :
60 : PRIVATE
61 :
62 : ! Global parameters in this module
63 :
64 : CHARACTER(LEN=*), PARAMETER, PRIVATE :: rng_record_format = "(A40,I2,3L2,ES25.16,18F20.1)"
65 : INTEGER, PARAMETER :: rng_record_length = 433
66 : INTEGER, PARAMETER :: rng_name_length = 40
67 :
68 : ! Distribution types:
69 :
70 : ! GAUSSIAN: Gaussian distribution with zero mean and unit variance
71 : ! UNIFORM: Uniform distribution [0,1] with 1/2 mean (default)
72 :
73 : INTEGER, PARAMETER :: GAUSSIAN = 1, &
74 : UNIFORM = 2
75 :
76 : REAL(KIND=dp), PARAMETER :: norm = 2.328306549295727688e-10_dp, &
77 : m1 = 4294967087.0_dp, &
78 : m2 = 4294944443.0_dp, &
79 : a12 = 1403580.0_dp, &
80 : a13n = 810728.0_dp, &
81 : a21 = 527612.0_dp, &
82 : a23n = 1370589.0_dp, &
83 : two17 = 131072.0_dp, & ! 2**17
84 : two53 = 9007199254740992.0_dp, & ! 2**53
85 : fact = 5.9604644775390625e-8_dp ! 1/2**24
86 :
87 : !&<
88 : ! The following are the transition matrices of the two MRG components
89 : ! (in matrix form), raised to the powers 1, 2**76, 2**127, and -1
90 :
91 : ! Transition matrix for the first component raised to the power 2**0
92 : REAL(KIND=dp), DIMENSION(3, 3), PARAMETER :: a1p0 = RESHAPE([ &
93 : 0.0_dp, 0.0_dp, -810728.0_dp, &
94 : 1.0_dp, 0.0_dp, 1403580.0_dp, &
95 : 0.0_dp, 1.0_dp, 0.0_dp &
96 : ], [3,3])
97 :
98 : ! Transition matrix for the second component raised to the power 2**0
99 : REAL(KIND=dp), DIMENSION(3, 3), PARAMETER :: a2p0 = RESHAPE([ &
100 : 0.0_dp, 0.0_dp, -1370589.0_dp, &
101 : 1.0_dp, 0.0_dp, 0.0_dp, &
102 : 0.0_dp, 1.0_dp, 527612.0_dp &
103 : ], [3,3])
104 :
105 : ! Transition matrix for the first component raised to the power 2**76
106 : REAL(KIND=dp), DIMENSION(3, 3), PARAMETER :: a1p76 = RESHAPE([ &
107 : 82758667.0_dp, 3672831523.0_dp, 3672091415.0_dp, &
108 : 1871391091.0_dp, 69195019.0_dp, 3528743235.0_dp, &
109 : 4127413238.0_dp, 1871391091.0_dp, 69195019.0_dp &
110 : ], [3,3])
111 :
112 : ! Transition matrix for the second component raised to the power 2**76
113 : REAL(KIND=dp), DIMENSION(3, 3), PARAMETER :: a2p76 = RESHAPE([ &
114 : 1511326704.0_dp, 4292754251.0_dp, 3859662829.0_dp, &
115 : 3759209742.0_dp, 1511326704.0_dp, 4292754251.0_dp, &
116 : 1610795712.0_dp, 3889917532.0_dp, 3708466080.0_dp &
117 : ], [3,3])
118 :
119 : ! Transition matrix for the first component raised to the power 2**127
120 : REAL(KIND=dp), DIMENSION(3, 3), PARAMETER :: a1p127 = RESHAPE([ &
121 : 2427906178.0_dp, 226153695.0_dp, 1988835001.0_dp, &
122 : 3580155704.0_dp, 1230515664.0_dp, 986791581.0_dp, &
123 : 949770784.0_dp, 3580155704.0_dp, 1230515664.0_dp &
124 : ], [3,3])
125 :
126 : ! Transition matrix for the second component raised to the power 2**127
127 : REAL(KIND=dp), DIMENSION(3, 3), PARAMETER :: a2p127 = RESHAPE([ &
128 : 1464411153.0_dp, 32183930.0_dp, 2824425944.0_dp, &
129 : 277697599.0_dp, 1464411153.0_dp, 32183930.0_dp, &
130 : 1610723613.0_dp, 1022607788.0_dp, 2093834863.0_dp &
131 : ], [3,3])
132 :
133 : ! Inverse of a1p0
134 : REAL(KIND=dp), DIMENSION(3, 3), PARAMETER :: inv_a1 = RESHAPE([ &
135 : 184888585.0_dp, 1.0_dp, 0.0_dp, &
136 : 0.0_dp, 0.0_dp, 1.0_dp, &
137 : 1945170933.0_dp, 0.0_dp, 0.0_dp &
138 : ], [3,3])
139 :
140 : ! Inverse of a2p0
141 : REAL(KIND=dp), DIMENSION(3, 3), PARAMETER :: inv_a2 = RESHAPE([ &
142 : 0.0_dp, 1.0_dp, 0.0_dp, &
143 : 360363334.0_dp, 0.0_dp, 1.0_dp, &
144 : 4225571728.0_dp, 0.0_dp, 0.0_dp &
145 : ], [3,3])
146 : !&>
147 :
148 : ! Data type definitions
149 :
150 : ! Information on a stream. The arrays bg, cg, and ig contain the current
151 : ! state of the stream, the starting state of the current substream, and the
152 : ! starting state of the stream. This stream generates antithetic variates
153 : ! if antithetic = .TRUE.. It also generates numbers with extended precision
154 : ! (53 bits, if machine follows IEEE 754 standard), if extended_precision =
155 : ! .TRUE., otherwise, numbers with 32 bits precision are generated.
156 :
157 : TYPE rng_stream_type
158 : PRIVATE
159 : ! the name could be an allocatable, but gfortran (even with 9.1) does not properly implement
160 : ! automatic deallocation of it and a `final`routine which would do it triggers an ICE in 7.4.1
161 : CHARACTER(LEN=rng_name_length) :: name = ""
162 : INTEGER :: distribution_type = UNIFORM
163 : ! ig: initial state, cg: current state, bg: initial state of the substream
164 : REAL(KIND=dp), DIMENSION(3, 2) :: bg = 0.0_dp, cg = 0.0_dp, ig = 0.0_dp
165 : LOGICAL :: antithetic = .FALSE., extended_precision = .FALSE.
166 : ! only used for distribution type GAUSSIAN
167 : REAL(KIND=dp) :: buffer = 0.0_dp
168 : LOGICAL :: buffer_filled = .FALSE.
169 :
170 : CONTAINS
171 : PROCEDURE, PASS(self) :: fill_1
172 : PROCEDURE, PASS(self) :: fill_2
173 : PROCEDURE, PASS(self) :: fill_3
174 : GENERIC, PUBLIC :: fill => fill_1, fill_2, fill_3
175 :
176 : PROCEDURE, PASS(self) :: next_int
177 : PROCEDURE, PASS(self) :: next_real
178 : GENERIC, PUBLIC :: next => next_int, next_real
179 :
180 : PROCEDURE, PASS(self), PUBLIC :: dump
181 : PROCEDURE, PASS(self), PUBLIC :: write
182 : PROCEDURE, PASS(self), PUBLIC :: advance
183 : PROCEDURE, PASS(self), PUBLIC :: set
184 : PROCEDURE, PASS(self), PUBLIC :: get
185 : PROCEDURE, PASS(self), PUBLIC :: reset
186 : PROCEDURE, PASS(self), PUBLIC :: reset_to_substream
187 : PROCEDURE, PASS(self), PUBLIC :: reset_to_next_substream
188 : PROCEDURE, PASS(self), PUBLIC :: shuffle
189 : END TYPE rng_stream_type
190 :
191 : INTERFACE rng_stream_type
192 : MODULE PROCEDURE :: rng_stream_constructor
193 : END INTERFACE
194 :
195 : TYPE rng_stream_p_type
196 : TYPE(rng_stream_type), POINTER :: stream => NULL()
197 : END TYPE rng_stream_p_type
198 :
199 : ! Public parameters
200 :
201 : PUBLIC :: rng_record_length, &
202 : rng_name_length, &
203 : GAUSSIAN, &
204 : UNIFORM
205 :
206 : ! Public data types
207 :
208 : PUBLIC :: rng_stream_p_type, &
209 : rng_stream_type
210 :
211 : ! Public subroutines
212 :
213 : PUBLIC :: check_rng, &
214 : next_rng_seed, &
215 : write_rng_matrices, &
216 : rng_stream_type_from_record
217 :
218 : CONTAINS
219 :
220 : ! **************************************************************************************************
221 : !> \brief Advance the state by n steps, i.e. jump n steps forward, if n > 0, or backward if n < 0.
222 : !> \param self ...
223 : !> \param e IF e > 0, let n = 2**e + c, IF e < 0, let n = -2**(-e) + c, IF e = 0, let n = c
224 : !> \param c ...
225 : !> \note The use of this method is discouraged
226 : ! **************************************************************************************************
227 112 : SUBROUTINE advance(self, e, c)
228 : CLASS(rng_stream_type), INTENT(INOUT) :: self
229 : INTEGER, INTENT(IN) :: e, c
230 :
231 : REAL(KIND=dp), DIMENSION(3, 2) :: x
232 : REAL(KIND=dp), DIMENSION(3, 3) :: u1, u2, v1, v2, w1, w2
233 :
234 112 : u1 = 0.0_dp
235 112 : u2 = 0.0_dp
236 112 : v1 = 0.0_dp
237 112 : v2 = 0.0_dp
238 112 : w1 = 0.0_dp
239 112 : w2 = 0.0_dp
240 :
241 112 : IF (e > 0) THEN
242 4 : CALL mat_two_pow_mod_m(a1p0, u1, m1, e)
243 4 : CALL mat_two_pow_mod_m(a2p0, u2, m2, e)
244 108 : ELSE IF (e < 0) THEN
245 3 : CALL mat_two_pow_mod_m(inv_a1, u1, m1, -e)
246 3 : CALL mat_two_pow_mod_m(inv_a2, u2, m2, -e)
247 : END IF
248 :
249 112 : IF (c >= 0) THEN
250 112 : CALL mat_pow_mod_m(a1p0, v1, m1, c)
251 112 : CALL mat_pow_mod_m(a2p0, v2, m2, c)
252 : ELSE
253 0 : CALL mat_pow_mod_m(inv_a1, v1, m1, -c)
254 0 : CALL mat_pow_mod_m(inv_a2, v2, m2, -c)
255 : END IF
256 :
257 112 : IF (e == 0) THEN
258 105 : w1 = v1
259 105 : w2 = v2
260 : ELSE
261 7 : CALL mat_mat_mod_m(u1, v1, w1, m1)
262 7 : CALL mat_mat_mod_m(u2, v2, w2, m2)
263 : END IF
264 :
265 112 : x = 0.0_dp
266 :
267 112 : CALL mat_vec_mod_m(w1, self%cg(:, 1), x(:, 1), m1)
268 112 : CALL mat_vec_mod_m(w2, self%cg(:, 2), x(:, 2), m2)
269 :
270 1008 : self%cg = x
271 112 : END SUBROUTINE advance
272 :
273 : ! **************************************************************************************************
274 : !> \brief ...
275 : !> \param output_unit ...
276 : !> \param ionode ...
277 : ! **************************************************************************************************
278 3 : SUBROUTINE check_rng(output_unit, ionode)
279 :
280 : ! Check the parallel (pseudo)random number generator (RNG).
281 :
282 : INTEGER, INTENT(IN) :: output_unit
283 : LOGICAL, INTENT(IN) :: ionode
284 :
285 : INTEGER :: i, sumi
286 : REAL(KIND=dp) :: sum, sum3
287 : REAL(KIND=dp), DIMENSION(3, 2) :: germe
288 : TYPE(rng_stream_type) :: cantor, g1, g2, g3, galois, laplace, &
289 : poisson
290 :
291 : ! -------------------------------------------------------------------------
292 : ! Test 1
293 :
294 : ! Create RNG test streams
295 :
296 3 : g1 = rng_stream_type("g1")
297 3 : g2 = rng_stream_type("g2", g1)
298 3 : g3 = rng_stream_type("g3", g2)
299 :
300 3 : IF (ionode) THEN
301 : WRITE (UNIT=output_unit, FMT="(/,T2,A)") &
302 2 : "RESULTS OF THE PSEUDO(RANDOM) NUMBER GENERATOR TEST RUNS", &
303 4 : "Initial states of the (pseudo)random number streams (test 1):"
304 2 : CALL g1%write(output_unit)
305 2 : CALL g2%write(output_unit)
306 2 : CALL g3%write(output_unit)
307 : END IF
308 :
309 3 : sum = g2%next() + g3%next()
310 :
311 3 : CALL g1%advance(5, 3)
312 3 : sum = sum + g1%next()
313 :
314 3 : CALL g1%reset()
315 108 : DO i = 1, 35
316 108 : CALL g1%advance(0, 1)
317 : END DO
318 3 : sum = sum + g1%next()
319 :
320 3 : CALL g1%reset()
321 :
322 3 : sumi = 0
323 108 : DO i = 1, 35
324 108 : sumi = sumi + g1%next(1, 10)
325 : END DO
326 3 : sum = sum + sumi/100.0_dp
327 :
328 3 : sum3 = 0.0_dp
329 303 : DO i = 1, 100
330 303 : sum3 = sum3 + g3%next()
331 : END DO
332 3 : sum = sum + sum3/10.0_dp
333 :
334 3 : CALL g3%reset()
335 18 : DO i = 1, 5
336 18 : sum = sum + g3%next()
337 : END DO
338 :
339 3 : CALL g3%reset()
340 15 : DO i = 1, 4
341 15 : CALL g3%reset_to_next_substream()
342 : END DO
343 18 : DO i = 1, 5
344 18 : sum = sum + g3%next()
345 : END DO
346 :
347 3 : CALL g3%reset_to_substream()
348 18 : DO i = 1, 5
349 18 : sum = sum + g3%next()
350 : END DO
351 :
352 3 : CALL g2%reset_to_next_substream()
353 3 : sum3 = 0.0_dp
354 300003 : DO i = 1, 100000
355 300003 : sum3 = sum3 + g2%next()
356 : END DO
357 3 : sum = sum + sum3/10000.0_dp
358 :
359 3 : CALL g3%set(antithetic=.TRUE.)
360 3 : sum3 = 0.0_dp
361 300003 : DO i = 1, 100000
362 300003 : sum3 = sum3 + g3%next()
363 : END DO
364 3 : sum = sum + sum3/10000.0_dp
365 :
366 3 : IF (ionode) THEN
367 : WRITE (UNIT=output_unit, FMT="(/,T2,A)") &
368 2 : "Final states of the (pseudo)random number streams (test 1):"
369 2 : CALL g1%write(output_unit)
370 2 : CALL g2%write(output_unit)
371 2 : CALL g3%write(output_unit)
372 : WRITE (UNIT=output_unit, FMT="(/,(T2,A))") &
373 2 : "This test routine should print for test 1 the number 25.342059"
374 : WRITE (UNIT=output_unit, FMT="(T2,A,F10.6)") &
375 2 : "The actual result of test 1 is ", sum
376 : END IF
377 :
378 : ! -------------------------------------------------------------------------
379 : ! Test 2
380 :
381 27 : germe(:, :) = 1
382 :
383 3 : poisson = rng_stream_type("Poisson", seed=germe)
384 3 : laplace = rng_stream_type("Laplace", poisson)
385 3 : galois = rng_stream_type("Galois", laplace)
386 3 : cantor = rng_stream_type("Cantor", galois)
387 :
388 3 : IF (ionode) THEN
389 : WRITE (UNIT=output_unit, FMT="(/,T2,A)") &
390 2 : "Initial states of the (pseudo)random number streams (test 2):"
391 2 : CALL poisson%write(output_unit)
392 2 : CALL laplace%write(output_unit)
393 2 : CALL galois%write(output_unit)
394 2 : CALL cantor%write(output_unit)
395 : END IF
396 :
397 3 : sum = sum + poisson%next() + laplace%next() + galois%next() + cantor%next()
398 :
399 3 : CALL galois%advance(-127, 0)
400 3 : sum = sum + galois%next()
401 :
402 3 : CALL galois%reset_to_next_substream()
403 3 : CALL galois%set(extended_precision=.TRUE.)
404 3 : sum3 = 0.0_dp
405 300003 : DO i = 1, 100000
406 300003 : sum3 = sum3 + galois%next()
407 : END DO
408 3 : sum = sum + sum3/10000.0_dp
409 :
410 3 : CALL galois%set(antithetic=.TRUE.)
411 3 : sum3 = 0.0_dp
412 300003 : DO i = 1, 100000
413 300003 : sum3 = sum3 + galois%next()
414 : END DO
415 3 : sum = sum + sum3/10000.0_dp
416 3 : CALL galois%set(antithetic=.FALSE.)
417 :
418 3 : CALL galois%set(extended_precision=.FALSE.)
419 3 : sum = sum + poisson%next() + laplace%next() + galois%next() + cantor%next()
420 :
421 3 : IF (ionode) THEN
422 : WRITE (UNIT=output_unit, FMT="(/,T2,A)") &
423 2 : "Final states of the (pseudo)random number streams (test 2):"
424 2 : CALL poisson%write(output_unit)
425 2 : CALL laplace%write(output_unit)
426 2 : CALL galois%write(output_unit)
427 2 : CALL cantor%write(output_unit)
428 : WRITE (UNIT=output_unit, FMT="(/,(T2,A))") &
429 2 : "This test routine should print for test 2 the number 39.697547"
430 : WRITE (UNIT=output_unit, FMT="(T2,A,F10.6)") &
431 2 : "The actual result of test 2 is ", sum
432 : END IF
433 :
434 564 : END SUBROUTINE check_rng
435 :
436 : ! **************************************************************************************************
437 : !> \brief Check that the seeds are legitimate values.
438 : !> \param seed ...
439 : ! **************************************************************************************************
440 111036 : SUBROUTINE check_seed(seed)
441 : REAL(KIND=dp), DIMENSION(3, 2), INTENT(IN) :: seed
442 :
443 : CHARACTER(LEN=*), PARAMETER :: fmtstr = "(A,I1,A,ES23.14,A,ES23.14)"
444 :
445 : CHARACTER(LEN=default_string_length) :: message
446 : INTEGER :: i
447 :
448 444144 : DO i = 1, 3
449 :
450 : ! Check condition: 0 <= seed(:,1) < m1
451 :
452 333108 : IF (seed(i, 1) < 0.0_dp) THEN
453 : WRITE (UNIT=message, FMT=fmtstr) &
454 0 : "seed(", i, ",1) = ", seed(i, 1), " < ", 0.0_dp
455 0 : CALL compress(message)
456 0 : CPABORT(message)
457 : END IF
458 333108 : IF (seed(i, 1) >= m1) THEN
459 : WRITE (UNIT=message, FMT=fmtstr) &
460 0 : "seed(", i, ",1) = ", seed(i, 1), " >= ", m1
461 0 : CALL compress(message)
462 0 : CPABORT(message)
463 : END IF
464 :
465 : ! Check condition: 0 <= seed(:,2) < m2
466 :
467 333108 : IF (seed(i, 2) < 0.0_dp) THEN
468 : WRITE (UNIT=message, FMT=fmtstr) &
469 0 : "seed(", i, ",2) = ", seed(i, 2), " < ", 0.0_dp
470 0 : CALL compress(message)
471 0 : CPABORT(message)
472 : END IF
473 444144 : IF (seed(i, 2) >= m2) THEN
474 : WRITE (UNIT=message, FMT=fmtstr) &
475 0 : "seed(", i, ",2) = ", seed(i, 2), " >= ", m2
476 0 : CALL compress(message)
477 0 : CPABORT(message)
478 : END IF
479 :
480 : END DO
481 :
482 : ! Check condition: first or second seed is 0
483 :
484 111036 : IF (ALL(seed(:, 1) < 1.0_dp)) THEN
485 0 : CPABORT("First seed = 0")
486 : END IF
487 :
488 111036 : IF (ALL(seed(:, 2) < 1.0_dp)) THEN
489 0 : CPABORT("Second seed = 0")
490 : END IF
491 :
492 111036 : END SUBROUTINE check_seed
493 :
494 : ! **************************************************************************************************
495 : !> \brief Create a new RNG stream.
496 : !> \param name ...
497 : !> \param last_rng_stream ...
498 : !> \param distribution_type ...
499 : !> \param seed ...
500 : !> \param antithetic ...
501 : !> \param extended_precision ...
502 : !> \return ...
503 : ! **************************************************************************************************
504 52865 : FUNCTION rng_stream_constructor(name, last_rng_stream, distribution_type, seed, antithetic, extended_precision) &
505 1321625 : RESULT(rng_stream)
506 :
507 : CHARACTER(LEN=*), INTENT(IN) :: name
508 : TYPE(rng_stream_type), INTENT(IN), OPTIONAL :: last_rng_stream
509 : INTEGER, INTENT(IN), OPTIONAL :: distribution_type
510 : REAL(KIND=dp), DIMENSION(3, 2), INTENT(IN), &
511 : OPTIONAL :: seed
512 : LOGICAL, INTENT(IN), OPTIONAL :: antithetic, extended_precision
513 : TYPE(rng_stream_type) :: rng_stream
514 :
515 52865 : IF (LEN_TRIM(name) > rng_name_length) THEN
516 0 : CPABORT("given random number generator name is too long")
517 : END IF
518 :
519 52865 : rng_stream%name = TRIM(name)
520 :
521 52865 : IF (PRESENT(seed)) THEN
522 52566 : CALL check_seed(seed)
523 473094 : rng_stream%ig = seed
524 299 : ELSE IF (PRESENT(last_rng_stream)) THEN
525 1188 : rng_stream%ig = next_rng_seed(last_rng_stream%ig)
526 : ELSE
527 1503 : rng_stream%ig = next_rng_seed()
528 : END IF
529 :
530 475785 : rng_stream%cg = rng_stream%ig
531 475785 : rng_stream%bg = rng_stream%ig
532 :
533 52865 : IF (PRESENT(distribution_type)) THEN
534 102525 : SELECT CASE (distribution_type)
535 : CASE (GAUSSIAN)
536 49687 : rng_stream%distribution_type = GAUSSIAN
537 : CASE (UNIFORM)
538 3151 : rng_stream%distribution_type = UNIFORM
539 : CASE DEFAULT
540 52838 : CPABORT("Invalid distribution type specified")
541 : END SELECT
542 27 : ELSE IF (PRESENT(last_rng_stream)) THEN
543 15 : rng_stream%distribution_type = last_rng_stream%distribution_type
544 : END IF
545 :
546 52865 : IF (PRESENT(antithetic)) THEN
547 0 : rng_stream%antithetic = antithetic
548 52865 : ELSE IF (PRESENT(last_rng_stream)) THEN
549 132 : rng_stream%antithetic = last_rng_stream%antithetic
550 : END IF
551 :
552 52865 : IF (PRESENT(extended_precision)) THEN
553 52680 : rng_stream%extended_precision = extended_precision
554 185 : ELSE IF (PRESENT(last_rng_stream)) THEN
555 35 : rng_stream%extended_precision = last_rng_stream%extended_precision
556 : END IF
557 1321625 : END FUNCTION rng_stream_constructor
558 :
559 : ! **************************************************************************************************
560 : !> \brief Create a RNG stream from a record given as an internal file (string).
561 : !> \param rng_record ...
562 : !> \return ...
563 : ! **************************************************************************************************
564 36400 : FUNCTION rng_stream_type_from_record(rng_record) RESULT(rng_stream)
565 : CHARACTER(LEN=rng_record_length), INTENT(IN) :: rng_record
566 : TYPE(rng_stream_type) :: rng_stream
567 :
568 : READ (UNIT=rng_record, FMT=rng_record_format) &
569 1300 : rng_stream%name, &
570 1300 : rng_stream%distribution_type, &
571 1300 : rng_stream%antithetic, &
572 1300 : rng_stream%extended_precision, &
573 1300 : rng_stream%buffer_filled, &
574 1300 : rng_stream%buffer, &
575 1300 : rng_stream%cg, &
576 1300 : rng_stream%bg, &
577 2600 : rng_stream%ig
578 35100 : END FUNCTION rng_stream_type_from_record
579 :
580 : ! **************************************************************************************************
581 : !> \brief Dump a RNG stream as a record given as an internal file (string).
582 : !> \param self ...
583 : !> \param rng_record ...
584 : ! **************************************************************************************************
585 16925 : SUBROUTINE dump(self, rng_record)
586 : CLASS(rng_stream_type), INTENT(IN) :: self
587 : CHARACTER(LEN=rng_record_length), INTENT(OUT) :: rng_record
588 :
589 16925 : rng_record = " "
590 : WRITE (UNIT=rng_record, FMT=rng_record_format) &
591 16925 : self%name, &
592 16925 : self%distribution_type, &
593 16925 : self%antithetic, &
594 16925 : self%extended_precision, &
595 16925 : self%buffer_filled, &
596 16925 : self%buffer, &
597 16925 : self%cg, &
598 16925 : self%bg, &
599 33850 : self%ig
600 16925 : END SUBROUTINE dump
601 :
602 : ! **************************************************************************************************
603 : !> \brief Get the components of a RNG stream.
604 : !> \param self ...
605 : !> \param name ...
606 : !> \param distribution_type ...
607 : !> \param bg ...
608 : !> \param cg ...
609 : !> \param ig ...
610 : !> \param antithetic ...
611 : !> \param extended_precision ...
612 : !> \param buffer ...
613 : !> \param buffer_filled ...
614 : !> \par History
615 : !> 2009-11-04 changed bg, cg and ig type from INTEGER, DIMENSION(3, 2)
616 : !> to REAL(KIND=dp), DIMENSION(3, 2) [lwalewski]
617 : !> 2009-11-09 getting the buffer and buffer_filled components
618 : !> added [lwalewski]
619 : ! **************************************************************************************************
620 16710 : SUBROUTINE get(self, name, distribution_type, bg, cg, ig, &
621 : antithetic, extended_precision, &
622 : buffer, buffer_filled)
623 :
624 : CLASS(rng_stream_type), INTENT(IN) :: self
625 : CHARACTER(LEN=rng_name_length), INTENT(OUT), OPTIONAL :: name
626 : INTEGER, INTENT(OUT), OPTIONAL :: distribution_type
627 : REAL(KIND=dp), DIMENSION(3, 2), INTENT(OUT), &
628 : OPTIONAL :: bg, cg, ig
629 : LOGICAL, INTENT(OUT), OPTIONAL :: antithetic, extended_precision
630 : REAL(KIND=dp), INTENT(OUT), OPTIONAL :: buffer
631 : LOGICAL, INTENT(OUT), OPTIONAL :: buffer_filled
632 :
633 16710 : IF (PRESENT(name)) name = self%name
634 16710 : IF (PRESENT(distribution_type)) THEN
635 0 : distribution_type = self%distribution_type
636 : END IF
637 126982 : IF (PRESENT(bg)) bg = self%bg
638 126982 : IF (PRESENT(cg)) cg = self%cg
639 150390 : IF (PRESENT(ig)) ig = self%ig
640 16710 : IF (PRESENT(antithetic)) antithetic = self%antithetic
641 16710 : IF (PRESENT(extended_precision)) THEN
642 0 : extended_precision = self%extended_precision
643 : END IF
644 16710 : IF (PRESENT(buffer)) buffer = self%buffer
645 16710 : IF (PRESENT(buffer_filled)) buffer_filled = self%buffer_filled
646 16710 : END SUBROUTINE get
647 :
648 : ! **************************************************************************************************
649 : !> \brief Returns c = MODULO(a*b,m)
650 : !> \param a ...
651 : !> \param b ...
652 : !> \param c ...
653 : !> \param m ...
654 : ! **************************************************************************************************
655 1064 : PURE SUBROUTINE mat_mat_mod_m(a, b, c, m)
656 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN) :: a, b
657 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(OUT) :: c
658 : REAL(KIND=dp), INTENT(IN) :: m
659 :
660 : INTEGER :: i
661 :
662 4256 : DO i = 1, 3
663 4256 : CALL mat_vec_mod_m(a, b(:, i), c(:, i), m)
664 : END DO
665 :
666 1064 : END SUBROUTINE mat_mat_mod_m
667 :
668 : ! **************************************************************************************************
669 : !> \brief Compute matrix b = MODULO(a**n,m)
670 : !> \param a ...
671 : !> \param b ...
672 : !> \param m ...
673 : !> \param n ...
674 : ! **************************************************************************************************
675 224 : PURE SUBROUTINE mat_pow_mod_m(a, b, m, n)
676 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN) :: a
677 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(OUT) :: b
678 : REAL(KIND=dp), INTENT(IN) :: m
679 : INTEGER, INTENT(IN) :: n
680 :
681 : INTEGER :: i
682 : REAL(KIND=dp), DIMENSION(3, 3) :: u, v, w
683 :
684 : ! Initialize: u = v = a; b = I
685 :
686 224 : w = a
687 :
688 224 : b(1, 1) = 1.0_dp
689 224 : b(2, 1) = 0.0_dp
690 224 : b(3, 1) = 0.0_dp
691 224 : b(1, 2) = 0.0_dp
692 224 : b(2, 2) = 1.0_dp
693 224 : b(3, 2) = 0.0_dp
694 224 : b(1, 3) = 0.0_dp
695 224 : b(2, 3) = 0.0_dp
696 224 : b(3, 3) = 1.0_dp
697 :
698 : ! Compute b = MODULO(a**n,m) using the binary decomposition of n
699 :
700 224 : i = n
701 :
702 16 : DO
703 240 : IF (MODULO(i, 2) /= 0) THEN
704 228 : u = w
705 228 : v = b
706 228 : CALL mat_mat_mod_m(u, v, b, m)
707 : END IF
708 240 : i = i/2
709 240 : IF (i == 0) EXIT
710 16 : u = w
711 16 : v = w
712 16 : CALL mat_mat_mod_m(u, v, w, m)
713 : END DO
714 224 : END SUBROUTINE mat_pow_mod_m
715 :
716 : ! **************************************************************************************************
717 : !> \brief Compute matrix b = MODULO(a**(2**e),m)
718 : !> \param a ...
719 : !> \param b ...
720 : !> \param m ...
721 : !> \param e ...
722 : ! **************************************************************************************************
723 14 : PURE SUBROUTINE mat_two_pow_mod_m(a, b, m, e)
724 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN) :: a
725 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(OUT) :: b
726 : REAL(KIND=dp), INTENT(IN) :: m
727 : INTEGER, INTENT(IN) :: e
728 :
729 : INTEGER :: i
730 : REAL(KIND=dp), DIMENSION(3, 3) :: u, v
731 :
732 14 : b = a
733 :
734 820 : DO i = 1, e
735 806 : u = b
736 806 : v = b
737 820 : CALL mat_mat_mod_m(u, v, b, m)
738 : END DO
739 :
740 14 : END SUBROUTINE mat_two_pow_mod_m
741 :
742 : ! **************************************************************************************************
743 : !> \brief Returns v = MODULO(a*s,m). Assumes that -m < s(i) < m.
744 : !> \param a ...
745 : !> \param s ...
746 : !> \param v ...
747 : !> \param m ...
748 : ! **************************************************************************************************
749 177940 : PURE SUBROUTINE mat_vec_mod_m(a, s, v, m)
750 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN) :: a
751 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: s
752 : REAL(KIND=dp), DIMENSION(3), INTENT(OUT) :: v
753 : REAL(KIND=dp), INTENT(IN) :: m
754 :
755 : INTEGER :: i, j
756 : REAL(KIND=dp) :: a1, a2, c
757 :
758 177940 : v = 0.0_dp
759 :
760 711760 : DO i = 1, 3
761 2313220 : DO j = 1, 3
762 1601460 : a2 = a(i, j)
763 1601460 : c = v(i)
764 1601460 : v(i) = a2*s(j) + c
765 1601460 : IF ((v(i) >= two53) .OR. (v(i) <= -two53)) THEN
766 1526473 : a1 = INT(a2/two17)
767 1526473 : a2 = a2 - a1*two17
768 1526473 : v(i) = a1*s(j)
769 1526473 : a1 = INT(v(i)/m)
770 1526473 : v(i) = v(i) - a1*m
771 1526473 : v(i) = v(i)*two17 + a2*s(j) + c
772 : END IF
773 1601460 : a1 = INT(v(i)/m)
774 1601460 : v(i) = v(i) - a1*m
775 2135280 : IF (v(i) < 0.0_dp) v(i) = v(i) + m
776 : END DO
777 : END DO
778 :
779 177940 : END SUBROUTINE mat_vec_mod_m
780 :
781 : ! **************************************************************************************************
782 : !> \brief Get the next integer random number between low and high from the stream
783 : !> \param self ...
784 : !> \param low ...
785 : !> \param high ...
786 : !> \return ...
787 : ! **************************************************************************************************
788 65742 : FUNCTION next_int(self, low, high) RESULT(u)
789 : CLASS(rng_stream_type), INTENT(INOUT) :: self
790 : INTEGER, INTENT(IN) :: low, high
791 : INTEGER :: u
792 :
793 : REAL(KIND=dp) :: r
794 :
795 65742 : CPASSERT(self%distribution_type == UNIFORM)
796 :
797 65742 : r = self%next_real()
798 65742 : u = low + INT(r*REAL(high - low + 1, dp))
799 65742 : END FUNCTION next_int
800 :
801 : ! **************************************************************************************************
802 : !> \brief Get the next real random number from the stream rng_stream.
803 : !> \param self ...
804 : !> \param variance variance of the Gaussian distribution (defaults to 1)
805 : !> \return ...
806 : ! **************************************************************************************************
807 46879621 : FUNCTION next_real(self, variance) RESULT(u)
808 : CLASS(rng_stream_type), INTENT(INOUT) :: self
809 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: variance
810 : REAL(KIND=dp) :: u
811 :
812 : REAL(KIND=dp) :: f, r, u1, u2, var
813 :
814 84166985 : SELECT CASE (self%distribution_type)
815 : CASE (GAUSSIAN)
816 37287364 : var = 1.0_dp
817 37287364 : IF (PRESENT(variance)) var = variance
818 : ! take the random number from the buffer, if the buffer is filled
819 37287364 : IF (self%buffer_filled) THEN
820 18633750 : u = SQRT(var)*self%buffer
821 18633750 : self%buffer_filled = .FALSE.
822 : ELSE
823 : DO
824 23744622 : IF (self%extended_precision) THEN
825 23744622 : u1 = 2.0_dp*rn53(self) - 1.0_dp
826 23744622 : u2 = 2.0_dp*rn53(self) - 1.0_dp
827 : ELSE
828 0 : u1 = 2.0_dp*rn32(self) - 1.0_dp
829 0 : u2 = 2.0_dp*rn32(self) - 1.0_dp
830 : END IF
831 23744622 : r = u1*u1 + u2*u2
832 23744622 : IF ((r > 0.0_dp) .AND. (r < 1.0_dp)) EXIT
833 : END DO
834 : ! Box-Muller transformation
835 18653614 : f = SQRT(-2.0_dp*LOG(r)/r)
836 18653614 : u = SQRT(var)*f*u1
837 : ! save the second random number for the next call
838 18653614 : self%buffer = f*u2
839 18653614 : self%buffer_filled = .TRUE.
840 : END IF
841 : CASE (UNIFORM)
842 46879621 : IF (self%extended_precision) THEN
843 8923613 : u = rn53(self)
844 : ELSE
845 668644 : u = rn32(self)
846 : END IF
847 : END SELECT
848 46879621 : END FUNCTION next_real
849 :
850 : ! **************************************************************************************************
851 : !> \brief Get the seed for the next RNG stream w.r.t. a given seed.
852 : !> \param seed If the optional argument seed is missing, then the default seed is returned.
853 : !> \return ...
854 : ! **************************************************************************************************
855 58813 : FUNCTION next_rng_seed(seed) RESULT(next_seed)
856 : REAL(KIND=dp), DIMENSION(3, 2), INTENT(IN), &
857 : OPTIONAL :: seed
858 : REAL(KIND=dp), DIMENSION(3, 2) :: next_seed
859 :
860 58813 : IF (PRESENT(seed)) THEN
861 58470 : CALL check_seed(seed)
862 58470 : CALL mat_vec_mod_m(a1p127, seed(:, 1), next_seed(:, 1), m1)
863 58470 : CALL mat_vec_mod_m(a2p127, seed(:, 2), next_seed(:, 2), m2)
864 : ELSE
865 3087 : next_seed = 12345.0_dp ! default seed
866 : END IF
867 :
868 58813 : END FUNCTION next_rng_seed
869 :
870 : ! **************************************************************************************************
871 : !> \brief Fill entity array with random numbers from the RNG stream rng_stream
872 : !> \param self ...
873 : !> \param array ...
874 : ! **************************************************************************************************
875 20043 : SUBROUTINE fill_1(self, array)
876 : CLASS(rng_stream_type), INTENT(INOUT) :: self
877 : REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: array
878 :
879 : INTEGER :: i
880 :
881 1420164 : DO i = 1, SIZE(array)
882 1420164 : array(i) = self%next()
883 : END DO
884 20043 : END SUBROUTINE fill_1
885 :
886 : ! **************************************************************************************************
887 : !> \brief ...
888 : !> \param self ...
889 : !> \param array ...
890 : ! **************************************************************************************************
891 0 : SUBROUTINE fill_2(self, array)
892 : CLASS(rng_stream_type), INTENT(INOUT) :: self
893 : REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: array
894 :
895 : INTEGER :: i, j
896 :
897 0 : DO j = 1, SIZE(array, 2)
898 0 : DO i = 1, SIZE(array, 1)
899 0 : array(i, j) = self%next()
900 : END DO
901 : END DO
902 0 : END SUBROUTINE fill_2
903 :
904 : ! **************************************************************************************************
905 : !> \brief ...
906 : !> \param self ...
907 : !> \param array ...
908 : ! **************************************************************************************************
909 0 : SUBROUTINE fill_3(self, array)
910 : CLASS(rng_stream_type), INTENT(INOUT) :: self
911 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(OUT) :: array
912 :
913 : INTEGER :: i, j, k
914 :
915 0 : DO k = 1, SIZE(array, 3)
916 0 : DO j = 1, SIZE(array, 2)
917 0 : DO i = 1, SIZE(array, 1)
918 0 : array(i, j, k) = self%next()
919 : END DO
920 : END DO
921 : END DO
922 0 : END SUBROUTINE fill_3
923 :
924 : ! **************************************************************************************************
925 : !> \brief Reset a random number stream to its initial state.
926 : !> \param self ...
927 : ! **************************************************************************************************
928 13 : SUBROUTINE reset(self)
929 : CLASS(rng_stream_type), INTENT(INOUT) :: self
930 :
931 117 : self%cg = self%ig
932 117 : self%bg = self%ig
933 13 : END SUBROUTINE reset
934 :
935 : ! **************************************************************************************************
936 : !> \brief Reset a random number stream to the beginning of its current substream.
937 : !> \param self ...
938 : ! **************************************************************************************************
939 3 : SUBROUTINE reset_to_substream(self)
940 : CLASS(rng_stream_type), INTENT(INOUT) :: self
941 :
942 27 : self%cg = self%bg
943 3 : END SUBROUTINE reset_to_substream
944 :
945 : ! **************************************************************************************************
946 : !> \brief Reset a random number stream to the beginning of its next substream.
947 : !> \param self ...
948 : ! **************************************************************************************************
949 28792 : SUBROUTINE reset_to_next_substream(self)
950 : CLASS(rng_stream_type), INTENT(INOUT) :: self
951 :
952 : REAL(KIND=dp), DIMENSION(3, 2) :: u
953 :
954 28792 : u = 0.0_dp
955 :
956 28792 : CALL mat_vec_mod_m(a1p76, self%bg(:, 1), u(:, 1), m1)
957 28792 : CALL mat_vec_mod_m(a2p76, self%bg(:, 2), u(:, 2), m2)
958 :
959 259128 : self%bg = u
960 259128 : self%cg = u
961 28792 : END SUBROUTINE reset_to_next_substream
962 :
963 : ! **************************************************************************************************
964 : !> \brief Generate the next random number with standard precision (32 bits)
965 : !> \param rng_stream ...
966 : !> \return ...
967 : ! **************************************************************************************************
968 113494358 : FUNCTION rn32(rng_stream) RESULT(u)
969 : TYPE(rng_stream_type) :: rng_stream
970 : REAL(KIND=dp) :: u
971 :
972 : INTEGER :: k
973 : REAL(KIND=dp) :: p1, p2
974 :
975 : ! Component 1
976 :
977 113494358 : p1 = a12*rng_stream%cg(2, 1) - a13n*rng_stream%cg(1, 1)
978 113494358 : k = INT(p1/m1)
979 113494358 : p1 = p1 - k*m1
980 113494358 : IF (p1 < 0.0_dp) p1 = p1 + m1
981 113494358 : rng_stream%cg(1, 1) = rng_stream%cg(2, 1)
982 113494358 : rng_stream%cg(2, 1) = rng_stream%cg(3, 1)
983 113494358 : rng_stream%cg(3, 1) = p1
984 :
985 : ! Component 2
986 :
987 113494358 : p2 = a21*rng_stream%cg(3, 2) - a23n*rng_stream%cg(1, 2)
988 113494358 : k = INT(p2/m2)
989 113494358 : p2 = p2 - k*m2
990 113494358 : IF (p2 < 0.0_dp) p2 = p2 + m2
991 113494358 : rng_stream%cg(1, 2) = rng_stream%cg(2, 2)
992 113494358 : rng_stream%cg(2, 2) = rng_stream%cg(3, 2)
993 113494358 : rng_stream%cg(3, 2) = p2
994 :
995 : ! Combination
996 :
997 113494358 : IF (p1 > p2) THEN
998 56723723 : u = (p1 - p2)*norm
999 : ELSE
1000 56770635 : u = (p1 - p2 + m1)*norm
1001 : END IF
1002 :
1003 113494358 : IF (rng_stream%antithetic) u = 1.0_dp - u
1004 :
1005 113494358 : END FUNCTION rn32
1006 :
1007 : ! **************************************************************************************************
1008 : !> \brief Generate the next random number with extended precision (53 bits)
1009 : !> \param rng_stream ...
1010 : !> \return ...
1011 : ! **************************************************************************************************
1012 56412857 : FUNCTION rn53(rng_stream) RESULT(u)
1013 : TYPE(rng_stream_type) :: rng_stream
1014 : REAL(KIND=dp) :: u
1015 :
1016 56412857 : u = rn32(rng_stream)
1017 :
1018 : ! Note: rn32 returns 1 - u in the antithetic case
1019 :
1020 56412857 : IF (rng_stream%antithetic) THEN
1021 300000 : u = u + (rn32(rng_stream) - 1.0_dp)*fact
1022 300000 : IF (u < 0.0_dp) u = u + 1.0_dp
1023 : ELSE
1024 56112857 : u = u + rn32(rng_stream)*fact
1025 56112857 : IF (u >= 1.0_dp) u = u - 1.0_dp
1026 : END IF
1027 56412857 : END FUNCTION rn53
1028 :
1029 : ! **************************************************************************************************
1030 : !> \brief Set the components of a RNG stream.
1031 : !> \param self ...
1032 : !> \param name ...
1033 : !> \param distribution_type ...
1034 : !> \param bg ...
1035 : !> \param cg ...
1036 : !> \param ig ...
1037 : !> \param seed ...
1038 : !> \param antithetic ...
1039 : !> \param extended_precision ...
1040 : !> \param buffer ...
1041 : !> \param buffer_filled ...
1042 : !> \par History
1043 : !> 2009-11-09 setting the buffer and buffer_filled components
1044 : !> added [lwalewski]
1045 : ! **************************************************************************************************
1046 13513 : SUBROUTINE set(self, name, distribution_type, bg, cg, ig, &
1047 : seed, antithetic, extended_precision, &
1048 : buffer, buffer_filled)
1049 :
1050 : ! NOTE: The manipulation of an active RNG stream is discouraged.
1051 :
1052 : CLASS(rng_stream_type), INTENT(INOUT) :: self
1053 : CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: name
1054 : INTEGER, INTENT(IN), OPTIONAL :: distribution_type
1055 : REAL(KIND=dp), DIMENSION(3, 2), INTENT(IN), &
1056 : OPTIONAL :: bg, cg, ig, seed
1057 : LOGICAL, INTENT(IN), OPTIONAL :: antithetic, extended_precision
1058 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: buffer
1059 : LOGICAL, INTENT(IN), OPTIONAL :: buffer_filled
1060 :
1061 0 : IF (PRESENT(name)) self%name = name
1062 13513 : IF (PRESENT(distribution_type)) THEN
1063 0 : self%distribution_type = distribution_type
1064 : END IF
1065 121497 : IF (PRESENT(bg)) self%bg = bg
1066 121497 : IF (PRESENT(cg)) self%cg = cg
1067 121497 : IF (PRESENT(ig)) self%ig = ig
1068 13513 : IF (PRESENT(seed)) THEN
1069 : ! Sets the initial seed of the stream to seed
1070 : ! NOTE: The use of this method is discouraged
1071 0 : CALL check_seed(seed)
1072 0 : self%ig = seed
1073 0 : self%cg = seed
1074 0 : self%bg = seed
1075 : END IF
1076 13513 : IF (PRESENT(antithetic)) self%antithetic = antithetic
1077 13513 : IF (PRESENT(extended_precision)) THEN
1078 6 : self%extended_precision = extended_precision
1079 : END IF
1080 13513 : IF (PRESENT(buffer)) self%buffer = buffer
1081 13513 : IF (PRESENT(buffer_filled)) self%buffer_filled = buffer_filled
1082 13513 : END SUBROUTINE set
1083 :
1084 : ! **************************************************************************************************
1085 : !> \brief Write the transformation matrices of the two MRG components (raised to the specified output)
1086 : !> \param output_unit ...
1087 : ! **************************************************************************************************
1088 1 : SUBROUTINE write_rng_matrices(output_unit)
1089 : INTEGER, INTENT(IN) :: output_unit
1090 :
1091 : CHARACTER(LEN=40) :: fmtstr
1092 : INTEGER :: i, j
1093 :
1094 : ! Print the transformation matrices for both components
1095 :
1096 : WRITE (UNIT=output_unit, FMT="(/,T2,A)") &
1097 : "TRANSFORMATION MATRICES FOR THE PARALLEL (PSEUDO)RANDOM NUMBER "// &
1098 1 : "GENERATOR"
1099 :
1100 1 : fmtstr = "(/,T4,A,/,/,(2X,3F14.1))"
1101 :
1102 : WRITE (UNIT=output_unit, FMT=fmtstr) &
1103 4 : "A1", ((a1p0(i, j), j=1, 3), i=1, 3)
1104 :
1105 : WRITE (UNIT=output_unit, FMT=fmtstr) &
1106 4 : "A2", ((a2p0(i, j), j=1, 3), i=1, 3)
1107 :
1108 : WRITE (UNIT=output_unit, FMT=fmtstr) &
1109 4 : "A1**(2**76)", ((a1p76(i, j), j=1, 3), i=1, 3)
1110 :
1111 : WRITE (UNIT=output_unit, FMT=fmtstr) &
1112 4 : "A2**(2**76)", ((a2p76(i, j), j=1, 3), i=1, 3)
1113 :
1114 : WRITE (UNIT=output_unit, FMT=fmtstr) &
1115 4 : "A1**(2**127)", ((a1p127(i, j), j=1, 3), i=1, 3)
1116 :
1117 : WRITE (UNIT=output_unit, FMT=fmtstr) &
1118 4 : "A2**(2**127)", ((a2p127(i, j), j=1, 3), i=1, 3)
1119 :
1120 1 : END SUBROUTINE write_rng_matrices
1121 :
1122 : ! **************************************************************************************************
1123 : !> \brief ...
1124 : !> \param self ...
1125 : !> \param output_unit ...
1126 : !> \param write_all if .TRUE., then print all stream informations (the default is .FALSE.).
1127 : ! **************************************************************************************************
1128 33 : SUBROUTINE write (self, output_unit, write_all)
1129 : CLASS(rng_stream_type), INTENT(IN) :: self
1130 : INTEGER, INTENT(IN) :: output_unit
1131 : LOGICAL, INTENT(IN), OPTIONAL :: write_all
1132 :
1133 : LOGICAL :: my_write_all
1134 :
1135 33 : my_write_all = .FALSE.
1136 :
1137 33 : IF (PRESENT(write_all)) THEN
1138 2 : my_write_all = write_all
1139 : END IF
1140 :
1141 : WRITE (UNIT=output_unit, FMT="(/,T2,A,/)") &
1142 33 : "Random number stream <"//TRIM(self%name)//">:"
1143 :
1144 36 : SELECT CASE (self%distribution_type)
1145 : CASE (GAUSSIAN)
1146 : WRITE (UNIT=output_unit, FMT="(T4,A)") &
1147 : "Distribution type: "// &
1148 3 : "Normal Gaussian distribution with zero mean"
1149 : CASE (UNIFORM)
1150 : WRITE (UNIT=output_unit, FMT="(T4,A)") &
1151 : "Distribution type: "// &
1152 33 : "Uniform distribution [0,1] with 1/2 mean"
1153 : END SELECT
1154 :
1155 33 : IF (self%antithetic) THEN
1156 2 : WRITE (UNIT=output_unit, FMT="(T4,A)") "Antithetic: yes"
1157 : ELSE
1158 31 : WRITE (UNIT=output_unit, FMT="(T4,A)") "Antithetic: no"
1159 : END IF
1160 :
1161 33 : IF (self%extended_precision) THEN
1162 5 : WRITE (UNIT=output_unit, FMT="(T4,A)") "Precision: 53 Bit"
1163 : ELSE
1164 28 : WRITE (UNIT=output_unit, FMT="(T4,A)") "Precision: 32 Bit"
1165 : END IF
1166 :
1167 33 : IF (my_write_all) THEN
1168 :
1169 : WRITE (UNIT=output_unit, FMT="(/,T4,A,/,/,(T4,A,3F20.1))") &
1170 2 : "Initial state of the stream:", &
1171 2 : "Component 1:", self%ig(:, 1), &
1172 4 : "Component 2:", self%ig(:, 2)
1173 :
1174 : WRITE (UNIT=output_unit, FMT="(/,T4,A,/,/,(T4,A,3F20.1))") &
1175 2 : "Initial state of the current substream:", &
1176 2 : "Component 1:", self%bg(:, 1), &
1177 4 : "Component 2:", self%bg(:, 2)
1178 :
1179 : END IF
1180 :
1181 : WRITE (UNIT=output_unit, FMT="(/,T4,A,/,/,(T4,A,3F20.1))") &
1182 33 : "Current state of the stream:", &
1183 33 : "Component 1:", self%cg(:, 1), &
1184 66 : "Component 2:", self%cg(:, 2)
1185 33 : END SUBROUTINE write
1186 :
1187 : ! **************************************************************************************************
1188 : !> \brief Shuffle an array of integers (using the Fisher-Yates shuffle)
1189 : !> \param self ...
1190 : !> \param arr the integer array to be shuffled
1191 : ! **************************************************************************************************
1192 76 : SUBROUTINE shuffle(self, arr)
1193 : CLASS(rng_stream_type), INTENT(INOUT) :: self
1194 : INTEGER, DIMENSION(:), INTENT(INOUT) :: arr
1195 :
1196 : INTEGER :: idxa, idxb, tmp
1197 :
1198 196 : DO idxa = UBOUND(arr, 1), LBOUND(arr, 1) + 1, -1
1199 44 : idxb = self%next(LBOUND(arr, 1), idxa)
1200 44 : tmp = arr(idxa)
1201 44 : arr(idxa) = arr(idxb)
1202 120 : arr(idxb) = tmp
1203 : END DO
1204 76 : END SUBROUTINE shuffle
1205 :
1206 0 : END MODULE parallel_rng_types
|
