(git:374b731)
Loading...
Searching...
No Matches
parallel_rng_types.F
Go to the documentation of this file.
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! **************************************************************************************************
51MODULE 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
218CONTAINS
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 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 u1 = 0.0_dp
235 u2 = 0.0_dp
236 v1 = 0.0_dp
237 v2 = 0.0_dp
238 w1 = 0.0_dp
239 w2 = 0.0_dp
240
241 IF (e > 0) THEN
242 CALL mat_two_pow_mod_m(a1p0, u1, m1, e)
243 CALL mat_two_pow_mod_m(a2p0, u2, m2, e)
244 ELSE IF (e < 0) THEN
245 CALL mat_two_pow_mod_m(inv_a1, u1, m1, -e)
246 CALL mat_two_pow_mod_m(inv_a2, u2, m2, -e)
247 END IF
248
249 IF (c >= 0) THEN
250 CALL mat_pow_mod_m(a1p0, v1, m1, c)
251 CALL mat_pow_mod_m(a2p0, v2, m2, c)
252 ELSE
253 CALL mat_pow_mod_m(inv_a1, v1, m1, -c)
254 CALL mat_pow_mod_m(inv_a2, v2, m2, -c)
255 END IF
256
257 IF (e == 0) THEN
258 w1 = v1
259 w2 = v2
260 ELSE
261 CALL mat_mat_mod_m(u1, v1, w1, m1)
262 CALL mat_mat_mod_m(u2, v2, w2, m2)
263 END IF
264
265 x = 0.0_dp
266
267 CALL mat_vec_mod_m(w1, self%cg(:, 1), x(:, 1), m1)
268 CALL mat_vec_mod_m(w2, self%cg(:, 2), x(:, 2), m2)
269
270 self%cg = x
271 END SUBROUTINE advance
272
273! **************************************************************************************************
274!> \brief ...
275!> \param output_unit ...
276!> \param ionode ...
277! **************************************************************************************************
278 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 g1 = rng_stream_type("g1")
297 g2 = rng_stream_type("g2", g1)
298 g3 = rng_stream_type("g3", g2)
299
300 IF (ionode) THEN
301 WRITE (unit=output_unit, fmt="(/,T2,A)") &
302 "RESULTS OF THE PSEUDO(RANDOM) NUMBER GENERATOR TEST RUNS", &
303 "Initial states of the (pseudo)random number streams (test 1):"
304 CALL g1%write(output_unit)
305 CALL g2%write(output_unit)
306 CALL g3%write(output_unit)
307 END IF
308
309 sum = g2%next() + g3%next()
310
311 CALL g1%advance(5, 3)
312 sum = sum + g1%next()
313
314 CALL g1%reset()
315 DO i = 1, 35
316 CALL g1%advance(0, 1)
317 END DO
318 sum = sum + g1%next()
319
320 CALL g1%reset()
321
322 sumi = 0
323 DO i = 1, 35
324 sumi = sumi + g1%next(1, 10)
325 END DO
326 sum = sum + sumi/100.0_dp
327
328 sum3 = 0.0_dp
329 DO i = 1, 100
330 sum3 = sum3 + g3%next()
331 END DO
332 sum = sum + sum3/10.0_dp
333
334 CALL g3%reset()
335 DO i = 1, 5
336 sum = sum + g3%next()
337 END DO
338
339 CALL g3%reset()
340 DO i = 1, 4
341 CALL g3%reset_to_next_substream()
342 END DO
343 DO i = 1, 5
344 sum = sum + g3%next()
345 END DO
346
347 CALL g3%reset_to_substream()
348 DO i = 1, 5
349 sum = sum + g3%next()
350 END DO
351
352 CALL g2%reset_to_next_substream()
353 sum3 = 0.0_dp
354 DO i = 1, 100000
355 sum3 = sum3 + g2%next()
356 END DO
357 sum = sum + sum3/10000.0_dp
358
359 CALL g3%set(antithetic=.true.)
360 sum3 = 0.0_dp
361 DO i = 1, 100000
362 sum3 = sum3 + g3%next()
363 END DO
364 sum = sum + sum3/10000.0_dp
365
366 IF (ionode) THEN
367 WRITE (unit=output_unit, fmt="(/,T2,A)") &
368 "Final states of the (pseudo)random number streams (test 1):"
369 CALL g1%write(output_unit)
370 CALL g2%write(output_unit)
371 CALL g3%write(output_unit)
372 WRITE (unit=output_unit, fmt="(/,(T2,A))") &
373 "This test routine should print for test 1 the number 25.342059"
374 WRITE (unit=output_unit, fmt="(T2,A,F10.6)") &
375 "The actual result of test 1 is ", sum
376 END IF
377
378 ! -------------------------------------------------------------------------
379 ! Test 2
380
381 germe(:, :) = 1
382
383 poisson = rng_stream_type("Poisson", seed=germe)
384 laplace = rng_stream_type("Laplace", poisson)
385 galois = rng_stream_type("Galois", laplace)
386 cantor = rng_stream_type("Cantor", galois)
387
388 IF (ionode) THEN
389 WRITE (unit=output_unit, fmt="(/,T2,A)") &
390 "Initial states of the (pseudo)random number streams (test 2):"
391 CALL poisson%write(output_unit)
392 CALL laplace%write(output_unit)
393 CALL galois%write(output_unit)
394 CALL cantor%write(output_unit)
395 END IF
396
397 sum = sum + poisson%next() + laplace%next() + galois%next() + cantor%next()
398
399 CALL galois%advance(-127, 0)
400 sum = sum + galois%next()
401
402 CALL galois%reset_to_next_substream()
403 CALL galois%set(extended_precision=.true.)
404 sum3 = 0.0_dp
405 DO i = 1, 100000
406 sum3 = sum3 + galois%next()
407 END DO
408 sum = sum + sum3/10000.0_dp
409
410 CALL galois%set(antithetic=.true.)
411 sum3 = 0.0_dp
412 DO i = 1, 100000
413 sum3 = sum3 + galois%next()
414 END DO
415 sum = sum + sum3/10000.0_dp
416 CALL galois%set(antithetic=.false.)
417
418 CALL galois%set(extended_precision=.false.)
419 sum = sum + poisson%next() + laplace%next() + galois%next() + cantor%next()
420
421 IF (ionode) THEN
422 WRITE (unit=output_unit, fmt="(/,T2,A)") &
423 "Final states of the (pseudo)random number streams (test 2):"
424 CALL poisson%write(output_unit)
425 CALL laplace%write(output_unit)
426 CALL galois%write(output_unit)
427 CALL cantor%write(output_unit)
428 WRITE (unit=output_unit, fmt="(/,(T2,A))") &
429 "This test routine should print for test 2 the number 39.697547"
430 WRITE (unit=output_unit, fmt="(T2,A,F10.6)") &
431 "The actual result of test 2 is ", sum
432 END IF
433
434 END SUBROUTINE check_rng
435
436! **************************************************************************************************
437!> \brief Check that the seeds are legitimate values.
438!> \param seed ...
439! **************************************************************************************************
440 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 DO i = 1, 3
449
450 ! Check condition: 0 <= seed(:,1) < m1
451
452 IF (seed(i, 1) < 0.0_dp) THEN
453 WRITE (unit=message, fmt=fmtstr) &
454 "seed(", i, ",1) = ", seed(i, 1), " < ", 0.0_dp
455 CALL compress(message)
456 cpabort(message)
457 END IF
458 IF (seed(i, 1) >= m1) THEN
459 WRITE (unit=message, fmt=fmtstr) &
460 "seed(", i, ",1) = ", seed(i, 1), " >= ", m1
461 CALL compress(message)
462 cpabort(message)
463 END IF
464
465 ! Check condition: 0 <= seed(:,2) < m2
466
467 IF (seed(i, 2) < 0.0_dp) THEN
468 WRITE (unit=message, fmt=fmtstr) &
469 "seed(", i, ",2) = ", seed(i, 2), " < ", 0.0_dp
470 CALL compress(message)
471 cpabort(message)
472 END IF
473 IF (seed(i, 2) >= m2) THEN
474 WRITE (unit=message, fmt=fmtstr) &
475 "seed(", i, ",2) = ", seed(i, 2), " >= ", m2
476 CALL compress(message)
477 cpabort(message)
478 END IF
479
480 END DO
481
482 ! Check condition: first or second seed is 0
483
484 IF (all(seed(:, 1) < 1.0_dp)) THEN
485 cpabort("First seed = 0")
486 END IF
487
488 IF (all(seed(:, 2) < 1.0_dp)) THEN
489 cpabort("Second seed = 0")
490 END IF
491
492 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 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 IF (len_trim(name) .GT. rng_name_length) &
516 cpabort("given random number generator name is too long")
517
518 rng_stream%name = trim(name)
519
520 IF (PRESENT(seed)) THEN
521 CALL check_seed(seed)
522 rng_stream%ig = seed
523 ELSE IF (PRESENT(last_rng_stream)) THEN
524 rng_stream%ig = next_rng_seed(last_rng_stream%ig)
525 ELSE
526 rng_stream%ig = next_rng_seed()
527 END IF
528
529 rng_stream%cg = rng_stream%ig
530 rng_stream%bg = rng_stream%ig
531
532 IF (PRESENT(distribution_type)) THEN
533 SELECT CASE (distribution_type)
534 CASE (gaussian)
535 rng_stream%distribution_type = gaussian
536 CASE (uniform)
537 rng_stream%distribution_type = uniform
538 CASE DEFAULT
539 cpabort("Invalid distribution type specified")
540 END SELECT
541 ELSE IF (PRESENT(last_rng_stream)) THEN
542 rng_stream%distribution_type = last_rng_stream%distribution_type
543 END IF
544
545 IF (PRESENT(antithetic)) THEN
546 rng_stream%antithetic = antithetic
547 ELSE IF (PRESENT(last_rng_stream)) THEN
548 rng_stream%antithetic = last_rng_stream%antithetic
549 END IF
550
551 IF (PRESENT(extended_precision)) THEN
552 rng_stream%extended_precision = extended_precision
553 ELSE IF (PRESENT(last_rng_stream)) THEN
554 rng_stream%extended_precision = last_rng_stream%extended_precision
555 END IF
556 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 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 rng_stream%name, &
569 rng_stream%distribution_type, &
570 rng_stream%antithetic, &
571 rng_stream%extended_precision, &
572 rng_stream%buffer_filled, &
573 rng_stream%buffer, &
574 rng_stream%cg, &
575 rng_stream%bg, &
576 rng_stream%ig
577 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 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 rng_record = " "
589 WRITE (unit=rng_record, fmt=rng_record_format) &
590 self%name, &
591 self%distribution_type, &
592 self%antithetic, &
593 self%extended_precision, &
594 self%buffer_filled, &
595 self%buffer, &
596 self%cg, &
597 self%bg, &
598 self%ig
599 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 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 IF (PRESENT(name)) name = self%name
633 IF (PRESENT(distribution_type)) &
634 distribution_type = self%distribution_type
635 IF (PRESENT(bg)) bg = self%bg
636 IF (PRESENT(cg)) cg = self%cg
637 IF (PRESENT(ig)) ig = self%ig
638 IF (PRESENT(antithetic)) antithetic = self%antithetic
639 IF (PRESENT(extended_precision)) &
640 extended_precision = self%extended_precision
641 IF (PRESENT(buffer)) buffer = self%buffer
642 IF (PRESENT(buffer_filled)) buffer_filled = self%buffer_filled
643 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 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 DO i = 1, 3
660 CALL mat_vec_mod_m(a, b(:, i), c(:, i), m)
661 END DO
662
663 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 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 w = a
684
685 b(1, 1) = 1.0_dp
686 b(2, 1) = 0.0_dp
687 b(3, 1) = 0.0_dp
688 b(1, 2) = 0.0_dp
689 b(2, 2) = 1.0_dp
690 b(3, 2) = 0.0_dp
691 b(1, 3) = 0.0_dp
692 b(2, 3) = 0.0_dp
693 b(3, 3) = 1.0_dp
694
695 ! Compute b = MODULO(a**n,m) using the binary decomposition of n
696
697 i = n
698
699 DO
700 IF (modulo(i, 2) /= 0) THEN
701 u = w
702 v = b
703 CALL mat_mat_mod_m(u, v, b, m)
704 END IF
705 i = i/2
706 IF (i == 0) EXIT
707 u = w
708 v = w
709 CALL mat_mat_mod_m(u, v, w, m)
710 END DO
711 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 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 b = a
730
731 DO i = 1, e
732 u = b
733 v = b
734 CALL mat_mat_mod_m(u, v, b, m)
735 END DO
736
737 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 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 v = 0.0_dp
756
757 DO i = 1, 3
758 DO j = 1, 3
759 a2 = a(i, j)
760 c = v(i)
761 v(i) = a2*s(j) + c
762 IF ((v(i) >= two53) .OR. (v(i) <= -two53)) THEN
763 a1 = int(a2/two17)
764 a2 = a2 - a1*two17
765 v(i) = a1*s(j)
766 a1 = int(v(i)/m)
767 v(i) = v(i) - a1*m
768 v(i) = v(i)*two17 + a2*s(j) + c
769 END IF
770 a1 = int(v(i)/m)
771 v(i) = v(i) - a1*m
772 IF (v(i) < 0.0_dp) v(i) = v(i) + m
773 END DO
774 END DO
775
776 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 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 cpassert(self%distribution_type == uniform)
793
794 r = self%next_real()
795 u = low + int(r*real(high - low + 1, dp))
796 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 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 SELECT CASE (self%distribution_type)
812 CASE (gaussian)
813 var = 1.0_dp
814 IF (PRESENT(variance)) var = variance
815 ! take the random number from the buffer, if the buffer is filled
816 IF (self%buffer_filled) THEN
817 u = sqrt(var)*self%buffer
818 self%buffer_filled = .false.
819 ELSE
820 DO
821 IF (self%extended_precision) THEN
822 u1 = 2.0_dp*rn53(self) - 1.0_dp
823 u2 = 2.0_dp*rn53(self) - 1.0_dp
824 ELSE
825 u1 = 2.0_dp*rn32(self) - 1.0_dp
826 u2 = 2.0_dp*rn32(self) - 1.0_dp
827 END IF
828 r = u1*u1 + u2*u2
829 IF ((r > 0.0_dp) .AND. (r < 1.0_dp)) EXIT
830 END DO
831 ! Box-Muller transformation
832 f = sqrt(-2.0_dp*log(r)/r)
833 u = sqrt(var)*f*u1
834 ! save the second random number for the next call
835 self%buffer = f*u2
836 self%buffer_filled = .true.
837 END IF
838 CASE (uniform)
839 IF (self%extended_precision) THEN
840 u = rn53(self)
841 ELSE
842 u = rn32(self)
843 END IF
844 END SELECT
845 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 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 IF (PRESENT(seed)) THEN
858 CALL check_seed(seed)
859 CALL mat_vec_mod_m(a1p127, seed(:, 1), next_seed(:, 1), m1)
860 CALL mat_vec_mod_m(a2p127, seed(:, 2), next_seed(:, 2), m2)
861 ELSE
862 next_seed = 12345.0_dp ! default seed
863 END IF
864
865 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 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 DO i = 1, SIZE(array)
879 array(i) = self%next()
880 END DO
881 END SUBROUTINE fill_1
882
883! **************************************************************************************************
884!> \brief ...
885!> \param self ...
886!> \param array ...
887! **************************************************************************************************
888 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 DO j = 1, SIZE(array, 2)
895 DO i = 1, SIZE(array, 1)
896 array(i, j) = self%next()
897 END DO
898 END DO
899 END SUBROUTINE fill_2
900
901! **************************************************************************************************
902!> \brief ...
903!> \param self ...
904!> \param array ...
905! **************************************************************************************************
906 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 DO k = 1, SIZE(array, 3)
913 DO j = 1, SIZE(array, 2)
914 DO i = 1, SIZE(array, 1)
915 array(i, j, k) = self%next()
916 END DO
917 END DO
918 END DO
919 END SUBROUTINE fill_3
920
921! **************************************************************************************************
922!> \brief Reset a random number stream to its initial state.
923!> \param self ...
924! **************************************************************************************************
925 SUBROUTINE reset(self)
926 CLASS(rng_stream_type), INTENT(INOUT) :: self
927
928 self%cg = self%ig
929 self%bg = self%ig
930 END SUBROUTINE reset
931
932! **************************************************************************************************
933!> \brief Reset a random number stream to the beginning of its current substream.
934!> \param self ...
935! **************************************************************************************************
936 SUBROUTINE reset_to_substream(self)
937 CLASS(rng_stream_type), INTENT(INOUT) :: self
938
939 self%cg = self%bg
940 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 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 u = 0.0_dp
952
953 CALL mat_vec_mod_m(a1p76, self%bg(:, 1), u(:, 1), m1)
954 CALL mat_vec_mod_m(a2p76, self%bg(:, 2), u(:, 2), m2)
955
956 self%bg = u
957 self%cg = u
958 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 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 p1 = a12*rng_stream%cg(2, 1) - a13n*rng_stream%cg(1, 1)
975 k = int(p1/m1)
976 p1 = p1 - k*m1
977 IF (p1 < 0.0_dp) p1 = p1 + m1
978 rng_stream%cg(1, 1) = rng_stream%cg(2, 1)
979 rng_stream%cg(2, 1) = rng_stream%cg(3, 1)
980 rng_stream%cg(3, 1) = p1
981
982 ! Component 2
983
984 p2 = a21*rng_stream%cg(3, 2) - a23n*rng_stream%cg(1, 2)
985 k = int(p2/m2)
986 p2 = p2 - k*m2
987 IF (p2 < 0.0_dp) p2 = p2 + m2
988 rng_stream%cg(1, 2) = rng_stream%cg(2, 2)
989 rng_stream%cg(2, 2) = rng_stream%cg(3, 2)
990 rng_stream%cg(3, 2) = p2
991
992 ! Combination
993
994 IF (p1 > p2) THEN
995 u = (p1 - p2)*norm
996 ELSE
997 u = (p1 - p2 + m1)*norm
998 END IF
999
1000 IF (rng_stream%antithetic) u = 1.0_dp - u
1001
1002 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 FUNCTION rn53(rng_stream) RESULT(u)
1010 TYPE(rng_stream_type) :: rng_stream
1011 REAL(kind=dp) :: u
1012
1013 u = rn32(rng_stream)
1014
1015 ! Note: rn32 returns 1 - u in the antithetic case
1016
1017 IF (rng_stream%antithetic) THEN
1018 u = u + (rn32(rng_stream) - 1.0_dp)*fact
1019 IF (u < 0.0_dp) u = u + 1.0_dp
1020 ELSE
1021 u = u + rn32(rng_stream)*fact
1022 IF (u >= 1.0_dp) u = u - 1.0_dp
1023 END IF
1024 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 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 IF (PRESENT(name)) self%name = name
1059 IF (PRESENT(distribution_type)) THEN
1060 self%distribution_type = distribution_type
1061 END IF
1062 IF (PRESENT(bg)) self%bg = bg
1063 IF (PRESENT(cg)) self%cg = cg
1064 IF (PRESENT(ig)) self%ig = ig
1065 IF (PRESENT(seed)) THEN
1066 ! Sets the initial seed of the stream to seed
1067 ! NOTE: The use of this method is discouraged
1068 CALL check_seed(seed)
1069 self%ig = seed
1070 self%cg = seed
1071 self%bg = seed
1072 END IF
1073 IF (PRESENT(antithetic)) self%antithetic = antithetic
1074 IF (PRESENT(extended_precision)) THEN
1075 self%extended_precision = extended_precision
1076 END IF
1077 IF (PRESENT(buffer)) self%buffer = buffer
1078 IF (PRESENT(buffer_filled)) self%buffer_filled = buffer_filled
1079 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 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 "GENERATOR"
1096
1097 fmtstr = "(/,T4,A,/,/,(2X,3F14.1))"
1098
1099 WRITE (unit=output_unit, fmt=fmtstr) &
1100 "A1", ((a1p0(i, j), j=1, 3), i=1, 3)
1101
1102 WRITE (unit=output_unit, fmt=fmtstr) &
1103 "A2", ((a2p0(i, j), j=1, 3), i=1, 3)
1104
1105 WRITE (unit=output_unit, fmt=fmtstr) &
1106 "A1**(2**76)", ((a1p76(i, j), j=1, 3), i=1, 3)
1107
1108 WRITE (unit=output_unit, fmt=fmtstr) &
1109 "A2**(2**76)", ((a2p76(i, j), j=1, 3), i=1, 3)
1110
1111 WRITE (unit=output_unit, fmt=fmtstr) &
1112 "A1**(2**127)", ((a1p127(i, j), j=1, 3), i=1, 3)
1113
1114 WRITE (unit=output_unit, fmt=fmtstr) &
1115 "A2**(2**127)", ((a2p127(i, j), j=1, 3), i=1, 3)
1116
1117 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 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 my_write_all = .false.
1133
1134 IF (PRESENT(write_all)) &
1135 my_write_all = write_all
1136
1137 WRITE (unit=output_unit, fmt="(/,T2,A,/)") &
1138 "Random number stream <"//trim(self%name)//">:"
1139
1140 SELECT CASE (self%distribution_type)
1141 CASE (gaussian)
1142 WRITE (unit=output_unit, fmt="(T4,A)") &
1143 "Distribution type: "// &
1144 "Normal Gaussian distribution with zero mean"
1145 CASE (uniform)
1146 WRITE (unit=output_unit, fmt="(T4,A)") &
1147 "Distribution type: "// &
1148 "Uniform distribution [0,1] with 1/2 mean"
1149 END SELECT
1150
1151 IF (self%antithetic) THEN
1152 WRITE (unit=output_unit, fmt="(T4,A)") "Antithetic: yes"
1153 ELSE
1154 WRITE (unit=output_unit, fmt="(T4,A)") "Antithetic: no"
1155 END IF
1156
1157 IF (self%extended_precision) THEN
1158 WRITE (unit=output_unit, fmt="(T4,A)") "Precision: 53 Bit"
1159 ELSE
1160 WRITE (unit=output_unit, fmt="(T4,A)") "Precision: 32 Bit"
1161 END IF
1162
1163 IF (my_write_all) THEN
1164
1165 WRITE (unit=output_unit, fmt="(/,T4,A,/,/,(T4,A,3F20.1))") &
1166 "Initial state of the stream:", &
1167 "Component 1:", self%ig(:, 1), &
1168 "Component 2:", self%ig(:, 2)
1169
1170 WRITE (unit=output_unit, fmt="(/,T4,A,/,/,(T4,A,3F20.1))") &
1171 "Initial state of the current substream:", &
1172 "Component 1:", self%bg(:, 1), &
1173 "Component 2:", self%bg(:, 2)
1174
1175 END IF
1176
1177 WRITE (unit=output_unit, fmt="(/,T4,A,/,/,(T4,A,3F20.1))") &
1178 "Current state of the stream:", &
1179 "Component 1:", self%cg(:, 1), &
1180 "Component 2:", self%cg(:, 2)
1181 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 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 DO idxa = ubound(arr, 1), lbound(arr, 1) + 1, -1
1195 idxb = self%next(lbound(arr, 1), idxa)
1196 tmp = arr(idxa)
1197 arr(idxa) = arr(idxb)
1198 arr(idxb) = tmp
1199 END DO
1200 END SUBROUTINE
1201
1202END MODULE parallel_rng_types
modulo
static GRID_HOST_DEVICE int modulo(int a, int m)
Equivalent of Fortran's MODULO, which always return a positive number. https://gcc....
Definition grid_common.h:117
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
kinds::default_string_length
integer, parameter, public default_string_length
Definition kinds.F:57
parallel_rng_types
Parallel (pseudo)random number generator (RNG) for multiple streams and substreams of random numbers.
Definition parallel_rng_types.F:51
parallel_rng_types::advance
subroutine advance(self, e, c)
Advance the state by n steps, i.e. jump n steps forward, if n > 0, or backward if n < 0.
Definition parallel_rng_types.F:228
parallel_rng_types::rng_stream_type_from_record
type(rng_stream_type) function, public rng_stream_type_from_record(rng_record)
Create a RNG stream from a record given as an internal file (string).
Definition parallel_rng_types.F:564
parallel_rng_types::next_rng_seed
real(kind=dp) function, dimension(3, 2), public next_rng_seed(seed)
Get the seed for the next RNG stream w.r.t. a given seed.
Definition parallel_rng_types.F:853
parallel_rng_types::rng_name_length
integer, parameter, public rng_name_length
Definition parallel_rng_types.F:66
parallel_rng_types::rng_record_length
integer, parameter, public rng_record_length
Definition parallel_rng_types.F:65
parallel_rng_types::uniform
integer, parameter, public uniform
Definition parallel_rng_types.F:73
parallel_rng_types::check_rng
subroutine, public check_rng(output_unit, ionode)
...
Definition parallel_rng_types.F:279
parallel_rng_types::write_rng_matrices
subroutine, public write_rng_matrices(output_unit)
Write the transformation matrices of the two MRG components (raised to the specified output)
Definition parallel_rng_types.F:1086
parallel_rng_types::gaussian
integer, parameter, public gaussian
Definition parallel_rng_types.F:73
string_utilities
Utilities for string manipulations.
Definition string_utilities.F:16
string_utilities::compress
subroutine, public compress(string, full)
Eliminate multiple space characters in a string. If full is .TRUE., then all spaces are eliminated.
Definition string_utilities.F:3190
parallel_rng_types::rng_stream_p_type
Definition parallel_rng_types.F:195
parallel_rng_types::rng_stream_type
Definition parallel_rng_types.F:157