Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2024 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 549 : END SUBROUTINE check_rng
435 :
436 : ! **************************************************************************************************
437 : !> \brief Check that the seeds are legitimate values.
438 : !> \param seed ...
439 : ! **************************************************************************************************
440 108606 : 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 434424 : DO i = 1, 3
449 :
450 : ! Check condition: 0 <= seed(:,1) < m1
451 :
452 325818 : 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 325818 : 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 325818 : 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 434424 : 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 108606 : IF (ALL(seed(:, 1) < 1.0_dp)) THEN
485 0 : CPABORT("First seed = 0")
486 : END IF
487 :
488 108606 : IF (ALL(seed(:, 2) < 1.0_dp)) THEN
489 0 : CPABORT("Second seed = 0")
490 : END IF
491 :
492 108606 : 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 50415 : FUNCTION rng_stream_constructor(name, last_rng_stream, distribution_type, seed, antithetic, extended_precision) &
505 : 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 50415 : IF (LEN_TRIM(name) .GT. rng_name_length) &
516 0 : CPABORT("given random number generator name is too long")
517 :
518 50415 : rng_stream%name = TRIM(name)
519 :
520 50415 : IF (PRESENT(seed)) THEN
521 50140 : CALL check_seed(seed)
522 451260 : rng_stream%ig = seed
523 275 : ELSE IF (PRESENT(last_rng_stream)) THEN
524 1188 : rng_stream%ig = next_rng_seed(last_rng_stream%ig)
525 : ELSE
526 1287 : rng_stream%ig = next_rng_seed()
527 : END IF
528 :
529 453735 : rng_stream%cg = rng_stream%ig
530 453735 : rng_stream%bg = rng_stream%ig
531 :
532 50415 : IF (PRESENT(distribution_type)) THEN
533 97965 : SELECT CASE (distribution_type)
534 : CASE (GAUSSIAN)
535 47577 : rng_stream%distribution_type = GAUSSIAN
536 : CASE (UNIFORM)
537 2811 : rng_stream%distribution_type = UNIFORM
538 : CASE DEFAULT
539 50388 : CPABORT("Invalid distribution type specified")
540 : END SELECT
541 27 : ELSE IF (PRESENT(last_rng_stream)) THEN
542 15 : rng_stream%distribution_type = last_rng_stream%distribution_type
543 : END IF
544 :
545 50415 : IF (PRESENT(antithetic)) THEN
546 0 : rng_stream%antithetic = antithetic
547 50415 : ELSE IF (PRESENT(last_rng_stream)) THEN
548 132 : rng_stream%antithetic = last_rng_stream%antithetic
549 : END IF
550 :
551 50415 : IF (PRESENT(extended_precision)) THEN
552 50252 : rng_stream%extended_precision = extended_precision
553 163 : ELSE IF (PRESENT(last_rng_stream)) THEN
554 35 : rng_stream%extended_precision = last_rng_stream%extended_precision
555 : END IF
556 1260375 : END FUNCTION rng_stream_constructor
557 :
558 : ! **************************************************************************************************
559 : !> \brief Create a RNG stream from a record given as an internal file (string).
560 : !> \param rng_record ...
561 : !> \return ...
562 : ! **************************************************************************************************
563 1300 : FUNCTION rng_stream_type_from_record(rng_record) RESULT(rng_stream)
564 : CHARACTER(LEN=rng_record_length), INTENT(IN) :: rng_record
565 : TYPE(rng_stream_type) :: rng_stream
566 :
567 : READ (UNIT=rng_record, FMT=rng_record_format) &
568 1300 : rng_stream%name, &
569 1300 : rng_stream%distribution_type, &
570 1300 : rng_stream%antithetic, &
571 1300 : rng_stream%extended_precision, &
572 1300 : rng_stream%buffer_filled, &
573 1300 : rng_stream%buffer, &
574 1300 : rng_stream%cg, &
575 1300 : rng_stream%bg, &
576 1300 : rng_stream%ig
577 35100 : END FUNCTION rng_stream_type_from_record
578 :
579 : ! **************************************************************************************************
580 : !> \brief Dump a RNG stream as a record given as an internal file (string).
581 : !> \param self ...
582 : !> \param rng_record ...
583 : ! **************************************************************************************************
584 16721 : SUBROUTINE dump(self, rng_record)
585 : CLASS(rng_stream_type), INTENT(IN) :: self
586 : CHARACTER(LEN=rng_record_length), INTENT(OUT) :: rng_record
587 :
588 16721 : rng_record = " "
589 : WRITE (UNIT=rng_record, FMT=rng_record_format) &
590 16721 : self%name, &
591 16721 : self%distribution_type, &
592 16721 : self%antithetic, &
593 16721 : self%extended_precision, &
594 16721 : self%buffer_filled, &
595 16721 : self%buffer, &
596 16721 : self%cg, &
597 16721 : self%bg, &
598 33442 : self%ig
599 16721 : END SUBROUTINE dump
600 :
601 : ! **************************************************************************************************
602 : !> \brief Get the components of a RNG stream.
603 : !> \param self ...
604 : !> \param name ...
605 : !> \param distribution_type ...
606 : !> \param bg ...
607 : !> \param cg ...
608 : !> \param ig ...
609 : !> \param antithetic ...
610 : !> \param extended_precision ...
611 : !> \param buffer ...
612 : !> \param buffer_filled ...
613 : !> \par History
614 : !> 2009-11-04 changed bg, cg and ig type from INTEGER, DIMENSION(3, 2)
615 : !> to REAL(KIND=dp), DIMENSION(3, 2) [lwalewski]
616 : !> 2009-11-09 getting the buffer and buffer_filled components
617 : !> added [lwalewski]
618 : ! **************************************************************************************************
619 16638 : SUBROUTINE get(self, name, distribution_type, bg, cg, ig, &
620 : antithetic, extended_precision, &
621 : buffer, buffer_filled)
622 :
623 : CLASS(rng_stream_type), INTENT(IN) :: self
624 : CHARACTER(LEN=rng_name_length), INTENT(OUT), OPTIONAL :: name
625 : INTEGER, INTENT(OUT), OPTIONAL :: distribution_type
626 : REAL(KIND=dp), DIMENSION(3, 2), INTENT(OUT), &
627 : OPTIONAL :: bg, cg, ig
628 : LOGICAL, INTENT(OUT), OPTIONAL :: antithetic, extended_precision
629 : REAL(KIND=dp), INTENT(OUT), OPTIONAL :: buffer
630 : LOGICAL, INTENT(OUT), OPTIONAL :: buffer_filled
631 :
632 16638 : IF (PRESENT(name)) name = self%name
633 16638 : IF (PRESENT(distribution_type)) &
634 0 : distribution_type = self%distribution_type
635 128878 : IF (PRESENT(bg)) bg = self%bg
636 128878 : IF (PRESENT(cg)) cg = self%cg
637 149742 : IF (PRESENT(ig)) ig = self%ig
638 16638 : IF (PRESENT(antithetic)) antithetic = self%antithetic
639 16638 : IF (PRESENT(extended_precision)) &
640 0 : extended_precision = self%extended_precision
641 16638 : IF (PRESENT(buffer)) buffer = self%buffer
642 16638 : IF (PRESENT(buffer_filled)) buffer_filled = self%buffer_filled
643 16638 : END SUBROUTINE get
644 :
645 : ! **************************************************************************************************
646 : !> \brief Returns c = MODULO(a*b,m)
647 : !> \param a ...
648 : !> \param b ...
649 : !> \param c ...
650 : !> \param m ...
651 : ! **************************************************************************************************
652 1064 : PURE SUBROUTINE mat_mat_mod_m(a, b, c, m)
653 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN) :: a, b
654 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(OUT) :: c
655 : REAL(KIND=dp), INTENT(IN) :: m
656 :
657 : INTEGER :: i
658 :
659 4256 : DO i = 1, 3
660 4256 : CALL mat_vec_mod_m(a, b(:, i), c(:, i), m)
661 : END DO
662 :
663 1064 : END SUBROUTINE mat_mat_mod_m
664 :
665 : ! **************************************************************************************************
666 : !> \brief Compute matrix b = MODULO(a**n,m)
667 : !> \param a ...
668 : !> \param b ...
669 : !> \param m ...
670 : !> \param n ...
671 : ! **************************************************************************************************
672 224 : PURE SUBROUTINE mat_pow_mod_m(a, b, m, n)
673 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN) :: a
674 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(OUT) :: b
675 : REAL(KIND=dp), INTENT(IN) :: m
676 : INTEGER, INTENT(IN) :: n
677 :
678 : INTEGER :: i
679 : REAL(KIND=dp), DIMENSION(3, 3) :: u, v, w
680 :
681 : ! Initialize: u = v = a; b = I
682 :
683 224 : w = a
684 :
685 224 : b(1, 1) = 1.0_dp
686 224 : b(2, 1) = 0.0_dp
687 224 : b(3, 1) = 0.0_dp
688 224 : b(1, 2) = 0.0_dp
689 224 : b(2, 2) = 1.0_dp
690 224 : b(3, 2) = 0.0_dp
691 224 : b(1, 3) = 0.0_dp
692 224 : b(2, 3) = 0.0_dp
693 224 : b(3, 3) = 1.0_dp
694 :
695 : ! Compute b = MODULO(a**n,m) using the binary decomposition of n
696 :
697 224 : i = n
698 :
699 16 : DO
700 240 : IF (MODULO(i, 2) /= 0) THEN
701 228 : u = w
702 228 : v = b
703 228 : CALL mat_mat_mod_m(u, v, b, m)
704 : END IF
705 240 : i = i/2
706 240 : IF (i == 0) EXIT
707 16 : u = w
708 16 : v = w
709 16 : CALL mat_mat_mod_m(u, v, w, m)
710 : END DO
711 224 : END SUBROUTINE mat_pow_mod_m
712 :
713 : ! **************************************************************************************************
714 : !> \brief Compute matrix b = MODULO(a**(2**e),m)
715 : !> \param a ...
716 : !> \param b ...
717 : !> \param m ...
718 : !> \param e ...
719 : ! **************************************************************************************************
720 14 : PURE SUBROUTINE mat_two_pow_mod_m(a, b, m, e)
721 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN) :: a
722 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(OUT) :: b
723 : REAL(KIND=dp), INTENT(IN) :: m
724 : INTEGER, INTENT(IN) :: e
725 :
726 : INTEGER :: i
727 : REAL(KIND=dp), DIMENSION(3, 3) :: u, v
728 :
729 14 : b = a
730 :
731 820 : DO i = 1, e
732 806 : u = b
733 806 : v = b
734 820 : CALL mat_mat_mod_m(u, v, b, m)
735 : END DO
736 :
737 14 : END SUBROUTINE mat_two_pow_mod_m
738 :
739 : ! **************************************************************************************************
740 : !> \brief Returns v = MODULO(a*s,m). Assumes that -m < s(i) < m.
741 : !> \param a ...
742 : !> \param s ...
743 : !> \param v ...
744 : !> \param m ...
745 : ! **************************************************************************************************
746 175824 : PURE SUBROUTINE mat_vec_mod_m(a, s, v, m)
747 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN) :: a
748 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: s
749 : REAL(KIND=dp), DIMENSION(3), INTENT(OUT) :: v
750 : REAL(KIND=dp), INTENT(IN) :: m
751 :
752 : INTEGER :: i, j
753 : REAL(KIND=dp) :: a1, a2, c
754 :
755 175824 : v = 0.0_dp
756 :
757 703296 : DO i = 1, 3
758 2285712 : DO j = 1, 3
759 1582416 : a2 = a(i, j)
760 1582416 : c = v(i)
761 1582416 : v(i) = a2*s(j) + c
762 1582416 : IF ((v(i) >= two53) .OR. (v(i) <= -two53)) THEN
763 1513563 : a1 = INT(a2/two17)
764 1513563 : a2 = a2 - a1*two17
765 1513563 : v(i) = a1*s(j)
766 1513563 : a1 = INT(v(i)/m)
767 1513563 : v(i) = v(i) - a1*m
768 1513563 : v(i) = v(i)*two17 + a2*s(j) + c
769 : END IF
770 1582416 : a1 = INT(v(i)/m)
771 1582416 : v(i) = v(i) - a1*m
772 2109888 : IF (v(i) < 0.0_dp) v(i) = v(i) + m
773 : END DO
774 : END DO
775 :
776 175824 : END SUBROUTINE mat_vec_mod_m
777 :
778 : ! **************************************************************************************************
779 : !> \brief Get the next integer random number between low and high from the stream
780 : !> \param self ...
781 : !> \param low ...
782 : !> \param high ...
783 : !> \return ...
784 : ! **************************************************************************************************
785 65742 : FUNCTION next_int(self, low, high) RESULT(u)
786 : CLASS(rng_stream_type), INTENT(INOUT) :: self
787 : INTEGER, INTENT(IN) :: low, high
788 : INTEGER :: u
789 :
790 : REAL(KIND=dp) :: r
791 :
792 65742 : CPASSERT(self%distribution_type == UNIFORM)
793 :
794 65742 : r = self%next_real()
795 65742 : u = low + INT(r*REAL(high - low + 1, dp))
796 65742 : END FUNCTION next_int
797 :
798 : ! **************************************************************************************************
799 : !> \brief Get the next real random number from the stream rng_stream.
800 : !> \param self ...
801 : !> \param variance variance of the Gaussian distribution (defaults to 1)
802 : !> \return ...
803 : ! **************************************************************************************************
804 46735443 : FUNCTION next_real(self, variance) RESULT(u)
805 : CLASS(rng_stream_type), INTENT(INOUT) :: self
806 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: variance
807 : REAL(KIND=dp) :: u
808 :
809 : REAL(KIND=dp) :: f, r, u1, u2, var
810 :
811 83750101 : SELECT CASE (self%distribution_type)
812 : CASE (GAUSSIAN)
813 37014658 : var = 1.0_dp
814 37014658 : IF (PRESENT(variance)) var = variance
815 : ! take the random number from the buffer, if the buffer is filled
816 37014658 : IF (self%buffer_filled) THEN
817 18497399 : u = SQRT(var)*self%buffer
818 18497399 : self%buffer_filled = .FALSE.
819 : ELSE
820 : DO
821 23571173 : IF (self%extended_precision) THEN
822 23571173 : u1 = 2.0_dp*rn53(self) - 1.0_dp
823 23571173 : u2 = 2.0_dp*rn53(self) - 1.0_dp
824 : ELSE
825 0 : u1 = 2.0_dp*rn32(self) - 1.0_dp
826 0 : u2 = 2.0_dp*rn32(self) - 1.0_dp
827 : END IF
828 23571173 : r = u1*u1 + u2*u2
829 23571173 : IF ((r > 0.0_dp) .AND. (r < 1.0_dp)) EXIT
830 : END DO
831 : ! Box-Muller transformation
832 18517259 : f = SQRT(-2.0_dp*LOG(r)/r)
833 18517259 : u = SQRT(var)*f*u1
834 : ! save the second random number for the next call
835 18517259 : self%buffer = f*u2
836 18517259 : self%buffer_filled = .TRUE.
837 : END IF
838 : CASE (UNIFORM)
839 46735443 : IF (self%extended_precision) THEN
840 9051329 : u = rn53(self)
841 : ELSE
842 669456 : u = rn32(self)
843 : END IF
844 : END SELECT
845 46735443 : END FUNCTION next_real
846 :
847 : ! **************************************************************************************************
848 : !> \brief Get the seed for the next RNG stream w.r.t. a given seed.
849 : !> \param seed If the optional argument seed is missing, then the default seed is returned.
850 : !> \return ...
851 : ! **************************************************************************************************
852 58781 : FUNCTION next_rng_seed(seed) RESULT(next_seed)
853 : REAL(KIND=dp), DIMENSION(3, 2), INTENT(IN), &
854 : OPTIONAL :: seed
855 : REAL(KIND=dp), DIMENSION(3, 2) :: next_seed
856 :
857 58781 : IF (PRESENT(seed)) THEN
858 58466 : CALL check_seed(seed)
859 58466 : CALL mat_vec_mod_m(a1p127, seed(:, 1), next_seed(:, 1), m1)
860 58466 : CALL mat_vec_mod_m(a2p127, seed(:, 2), next_seed(:, 2), m2)
861 : ELSE
862 2835 : next_seed = 12345.0_dp ! default seed
863 : END IF
864 :
865 58781 : END FUNCTION next_rng_seed
866 :
867 : ! **************************************************************************************************
868 : !> \brief Fill entity array with random numbers from the RNG stream rng_stream
869 : !> \param self ...
870 : !> \param array ...
871 : ! **************************************************************************************************
872 18825 : SUBROUTINE fill_1(self, array)
873 : CLASS(rng_stream_type), INTENT(INOUT) :: self
874 : REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: array
875 :
876 : INTEGER :: i
877 :
878 1546662 : DO i = 1, SIZE(array)
879 1546662 : array(i) = self%next()
880 : END DO
881 18825 : END SUBROUTINE fill_1
882 :
883 : ! **************************************************************************************************
884 : !> \brief ...
885 : !> \param self ...
886 : !> \param array ...
887 : ! **************************************************************************************************
888 0 : SUBROUTINE fill_2(self, array)
889 : CLASS(rng_stream_type), INTENT(INOUT) :: self
890 : REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: array
891 :
892 : INTEGER :: i, j
893 :
894 0 : DO j = 1, SIZE(array, 2)
895 0 : DO i = 1, SIZE(array, 1)
896 0 : array(i, j) = self%next()
897 : END DO
898 : END DO
899 0 : END SUBROUTINE fill_2
900 :
901 : ! **************************************************************************************************
902 : !> \brief ...
903 : !> \param self ...
904 : !> \param array ...
905 : ! **************************************************************************************************
906 0 : SUBROUTINE fill_3(self, array)
907 : CLASS(rng_stream_type), INTENT(INOUT) :: self
908 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(OUT) :: array
909 :
910 : INTEGER :: i, j, k
911 :
912 0 : DO k = 1, SIZE(array, 3)
913 0 : DO j = 1, SIZE(array, 2)
914 0 : DO i = 1, SIZE(array, 1)
915 0 : array(i, j, k) = self%next()
916 : END DO
917 : END DO
918 : END DO
919 0 : END SUBROUTINE fill_3
920 :
921 : ! **************************************************************************************************
922 : !> \brief Reset a random number stream to its initial state.
923 : !> \param self ...
924 : ! **************************************************************************************************
925 13 : SUBROUTINE reset(self)
926 : CLASS(rng_stream_type), INTENT(INOUT) :: self
927 :
928 117 : self%cg = self%ig
929 117 : self%bg = self%ig
930 13 : END SUBROUTINE reset
931 :
932 : ! **************************************************************************************************
933 : !> \brief Reset a random number stream to the beginning of its current substream.
934 : !> \param self ...
935 : ! **************************************************************************************************
936 3 : SUBROUTINE reset_to_substream(self)
937 : CLASS(rng_stream_type), INTENT(INOUT) :: self
938 :
939 27 : self%cg = self%bg
940 3 : END SUBROUTINE reset_to_substream
941 :
942 : ! **************************************************************************************************
943 : !> \brief Reset a random number stream to the beginning of its next substream.
944 : !> \param self ...
945 : ! **************************************************************************************************
946 27738 : SUBROUTINE reset_to_next_substream(self)
947 : CLASS(rng_stream_type), INTENT(INOUT) :: self
948 :
949 : REAL(KIND=dp), DIMENSION(3, 2) :: u
950 :
951 27738 : u = 0.0_dp
952 :
953 27738 : CALL mat_vec_mod_m(a1p76, self%bg(:, 1), u(:, 1), m1)
954 27738 : CALL mat_vec_mod_m(a2p76, self%bg(:, 2), u(:, 2), m2)
955 :
956 249642 : self%bg = u
957 249642 : self%cg = u
958 27738 : END SUBROUTINE reset_to_next_substream
959 :
960 : ! **************************************************************************************************
961 : !> \brief Generate the next random number with standard precision (32 bits)
962 : !> \param rng_stream ...
963 : !> \return ...
964 : ! **************************************************************************************************
965 113056806 : FUNCTION rn32(rng_stream) RESULT(u)
966 : TYPE(rng_stream_type) :: rng_stream
967 : REAL(KIND=dp) :: u
968 :
969 : INTEGER :: k
970 : REAL(KIND=dp) :: p1, p2
971 :
972 : ! Component 1
973 :
974 113056806 : p1 = a12*rng_stream%cg(2, 1) - a13n*rng_stream%cg(1, 1)
975 113056806 : k = INT(p1/m1)
976 113056806 : p1 = p1 - k*m1
977 113056806 : IF (p1 < 0.0_dp) p1 = p1 + m1
978 113056806 : rng_stream%cg(1, 1) = rng_stream%cg(2, 1)
979 113056806 : rng_stream%cg(2, 1) = rng_stream%cg(3, 1)
980 113056806 : rng_stream%cg(3, 1) = p1
981 :
982 : ! Component 2
983 :
984 113056806 : p2 = a21*rng_stream%cg(3, 2) - a23n*rng_stream%cg(1, 2)
985 113056806 : k = INT(p2/m2)
986 113056806 : p2 = p2 - k*m2
987 113056806 : IF (p2 < 0.0_dp) p2 = p2 + m2
988 113056806 : rng_stream%cg(1, 2) = rng_stream%cg(2, 2)
989 113056806 : rng_stream%cg(2, 2) = rng_stream%cg(3, 2)
990 113056806 : rng_stream%cg(3, 2) = p2
991 :
992 : ! Combination
993 :
994 113056806 : IF (p1 > p2) THEN
995 56505777 : u = (p1 - p2)*norm
996 : ELSE
997 56551029 : u = (p1 - p2 + m1)*norm
998 : END IF
999 :
1000 113056806 : IF (rng_stream%antithetic) u = 1.0_dp - u
1001 :
1002 113056806 : END FUNCTION rn32
1003 :
1004 : ! **************************************************************************************************
1005 : !> \brief Generate the next random number with extended precision (53 bits)
1006 : !> \param rng_stream ...
1007 : !> \return ...
1008 : ! **************************************************************************************************
1009 56193675 : FUNCTION rn53(rng_stream) RESULT(u)
1010 : TYPE(rng_stream_type) :: rng_stream
1011 : REAL(KIND=dp) :: u
1012 :
1013 56193675 : u = rn32(rng_stream)
1014 :
1015 : ! Note: rn32 returns 1 - u in the antithetic case
1016 :
1017 56193675 : IF (rng_stream%antithetic) THEN
1018 300000 : u = u + (rn32(rng_stream) - 1.0_dp)*fact
1019 300000 : IF (u < 0.0_dp) u = u + 1.0_dp
1020 : ELSE
1021 55893675 : u = u + rn32(rng_stream)*fact
1022 55893675 : IF (u >= 1.0_dp) u = u - 1.0_dp
1023 : END IF
1024 56193675 : END FUNCTION rn53
1025 :
1026 : ! **************************************************************************************************
1027 : !> \brief Set the components of a RNG stream.
1028 : !> \param self ...
1029 : !> \param name ...
1030 : !> \param distribution_type ...
1031 : !> \param bg ...
1032 : !> \param cg ...
1033 : !> \param ig ...
1034 : !> \param seed ...
1035 : !> \param antithetic ...
1036 : !> \param extended_precision ...
1037 : !> \param buffer ...
1038 : !> \param buffer_filled ...
1039 : !> \par History
1040 : !> 2009-11-09 setting the buffer and buffer_filled components
1041 : !> added [lwalewski]
1042 : ! **************************************************************************************************
1043 13759 : SUBROUTINE set(self, name, distribution_type, bg, cg, ig, &
1044 : seed, antithetic, extended_precision, &
1045 : buffer, buffer_filled)
1046 :
1047 : ! NOTE: The manipulation of an active RNG stream is discouraged.
1048 :
1049 : CLASS(rng_stream_type), INTENT(INOUT) :: self
1050 : CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: name
1051 : INTEGER, INTENT(IN), OPTIONAL :: distribution_type
1052 : REAL(KIND=dp), DIMENSION(3, 2), INTENT(IN), &
1053 : OPTIONAL :: bg, cg, ig, seed
1054 : LOGICAL, INTENT(IN), OPTIONAL :: antithetic, extended_precision
1055 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: buffer
1056 : LOGICAL, INTENT(IN), OPTIONAL :: buffer_filled
1057 :
1058 0 : IF (PRESENT(name)) self%name = name
1059 13759 : IF (PRESENT(distribution_type)) THEN
1060 0 : self%distribution_type = distribution_type
1061 : END IF
1062 123711 : IF (PRESENT(bg)) self%bg = bg
1063 123711 : IF (PRESENT(cg)) self%cg = cg
1064 123711 : IF (PRESENT(ig)) self%ig = ig
1065 13759 : IF (PRESENT(seed)) THEN
1066 : ! Sets the initial seed of the stream to seed
1067 : ! NOTE: The use of this method is discouraged
1068 0 : CALL check_seed(seed)
1069 0 : self%ig = seed
1070 0 : self%cg = seed
1071 0 : self%bg = seed
1072 : END IF
1073 13759 : IF (PRESENT(antithetic)) self%antithetic = antithetic
1074 13759 : IF (PRESENT(extended_precision)) THEN
1075 6 : self%extended_precision = extended_precision
1076 : END IF
1077 13759 : IF (PRESENT(buffer)) self%buffer = buffer
1078 13759 : IF (PRESENT(buffer_filled)) self%buffer_filled = buffer_filled
1079 13759 : END SUBROUTINE set
1080 :
1081 : ! **************************************************************************************************
1082 : !> \brief Write the transformation matrices of the two MRG components (raised to the specified output)
1083 : !> \param output_unit ...
1084 : ! **************************************************************************************************
1085 1 : SUBROUTINE write_rng_matrices(output_unit)
1086 : INTEGER, INTENT(IN) :: output_unit
1087 :
1088 : CHARACTER(LEN=40) :: fmtstr
1089 : INTEGER :: i, j
1090 :
1091 : ! Print the transformation matrices for both components
1092 :
1093 : WRITE (UNIT=output_unit, FMT="(/,T2,A)") &
1094 : "TRANSFORMATION MATRICES FOR THE PARALLEL (PSEUDO)RANDOM NUMBER "// &
1095 1 : "GENERATOR"
1096 :
1097 1 : fmtstr = "(/,T4,A,/,/,(2X,3F14.1))"
1098 :
1099 : WRITE (UNIT=output_unit, FMT=fmtstr) &
1100 4 : "A1", ((a1p0(i, j), j=1, 3), i=1, 3)
1101 :
1102 : WRITE (UNIT=output_unit, FMT=fmtstr) &
1103 4 : "A2", ((a2p0(i, j), j=1, 3), i=1, 3)
1104 :
1105 : WRITE (UNIT=output_unit, FMT=fmtstr) &
1106 4 : "A1**(2**76)", ((a1p76(i, j), j=1, 3), i=1, 3)
1107 :
1108 : WRITE (UNIT=output_unit, FMT=fmtstr) &
1109 4 : "A2**(2**76)", ((a2p76(i, j), j=1, 3), i=1, 3)
1110 :
1111 : WRITE (UNIT=output_unit, FMT=fmtstr) &
1112 4 : "A1**(2**127)", ((a1p127(i, j), j=1, 3), i=1, 3)
1113 :
1114 : WRITE (UNIT=output_unit, FMT=fmtstr) &
1115 4 : "A2**(2**127)", ((a2p127(i, j), j=1, 3), i=1, 3)
1116 :
1117 1 : END SUBROUTINE write_rng_matrices
1118 :
1119 : ! **************************************************************************************************
1120 : !> \brief ...
1121 : !> \param self ...
1122 : !> \param output_unit ...
1123 : !> \param write_all if .TRUE., then print all stream informations (the default is .FALSE.).
1124 : ! **************************************************************************************************
1125 33 : SUBROUTINE write (self, output_unit, write_all)
1126 : CLASS(rng_stream_type), INTENT(IN) :: self
1127 : INTEGER, INTENT(IN) :: output_unit
1128 : LOGICAL, INTENT(IN), OPTIONAL :: write_all
1129 :
1130 : LOGICAL :: my_write_all
1131 :
1132 33 : my_write_all = .FALSE.
1133 :
1134 33 : IF (PRESENT(write_all)) &
1135 2 : my_write_all = write_all
1136 :
1137 : WRITE (UNIT=output_unit, FMT="(/,T2,A,/)") &
1138 33 : "Random number stream <"//TRIM(self%name)//">:"
1139 :
1140 36 : SELECT CASE (self%distribution_type)
1141 : CASE (GAUSSIAN)
1142 : WRITE (UNIT=output_unit, FMT="(T4,A)") &
1143 : "Distribution type: "// &
1144 3 : "Normal Gaussian distribution with zero mean"
1145 : CASE (UNIFORM)
1146 : WRITE (UNIT=output_unit, FMT="(T4,A)") &
1147 : "Distribution type: "// &
1148 33 : "Uniform distribution [0,1] with 1/2 mean"
1149 : END SELECT
1150 :
1151 33 : IF (self%antithetic) THEN
1152 2 : WRITE (UNIT=output_unit, FMT="(T4,A)") "Antithetic: yes"
1153 : ELSE
1154 31 : WRITE (UNIT=output_unit, FMT="(T4,A)") "Antithetic: no"
1155 : END IF
1156 :
1157 33 : IF (self%extended_precision) THEN
1158 5 : WRITE (UNIT=output_unit, FMT="(T4,A)") "Precision: 53 Bit"
1159 : ELSE
1160 28 : WRITE (UNIT=output_unit, FMT="(T4,A)") "Precision: 32 Bit"
1161 : END IF
1162 :
1163 33 : IF (my_write_all) THEN
1164 :
1165 : WRITE (UNIT=output_unit, FMT="(/,T4,A,/,/,(T4,A,3F20.1))") &
1166 2 : "Initial state of the stream:", &
1167 2 : "Component 1:", self%ig(:, 1), &
1168 4 : "Component 2:", self%ig(:, 2)
1169 :
1170 : WRITE (UNIT=output_unit, FMT="(/,T4,A,/,/,(T4,A,3F20.1))") &
1171 2 : "Initial state of the current substream:", &
1172 2 : "Component 1:", self%bg(:, 1), &
1173 4 : "Component 2:", self%bg(:, 2)
1174 :
1175 : END IF
1176 :
1177 : WRITE (UNIT=output_unit, FMT="(/,T4,A,/,/,(T4,A,3F20.1))") &
1178 33 : "Current state of the stream:", &
1179 33 : "Component 1:", self%cg(:, 1), &
1180 66 : "Component 2:", self%cg(:, 2)
1181 33 : END SUBROUTINE write
1182 :
1183 : ! **************************************************************************************************
1184 : !> \brief Shuffle an array of integers (using the Fisher-Yates shuffle)
1185 : !> \param self ...
1186 : !> \param arr the integer array to be shuffled
1187 : ! **************************************************************************************************
1188 54 : SUBROUTINE shuffle(self, arr)
1189 : CLASS(rng_stream_type), INTENT(INOUT) :: self
1190 : INTEGER, DIMENSION(:), INTENT(INOUT) :: arr
1191 :
1192 : INTEGER :: idxa, idxb, tmp
1193 :
1194 152 : DO idxa = UBOUND(arr, 1), LBOUND(arr, 1) + 1, -1
1195 44 : idxb = self%next(LBOUND(arr, 1), idxa)
1196 44 : tmp = arr(idxa)
1197 44 : arr(idxa) = arr(idxb)
1198 98 : arr(idxb) = tmp
1199 : END DO
1200 54 : END SUBROUTINE
1201 :
1202 0 : END MODULE parallel_rng_types
|
