(git:71c3ab0)
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-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! **************************************************************************************************
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) > rng_name_length) THEN
516 cpabort("given random number generator name is too long")
517 END IF
518
519 rng_stream%name = trim(name)
520
521 IF (PRESENT(seed)) THEN
522 CALL check_seed(seed)
523 rng_stream%ig = seed
524 ELSE IF (PRESENT(last_rng_stream)) THEN
525 rng_stream%ig = next_rng_seed(last_rng_stream%ig)
526 ELSE
527 rng_stream%ig = next_rng_seed()
528 END IF
529
530 rng_stream%cg = rng_stream%ig
531 rng_stream%bg = rng_stream%ig
532
533 IF (PRESENT(distribution_type)) THEN
534 SELECT CASE (distribution_type)
535 CASE (gaussian)
536 rng_stream%distribution_type = gaussian
537 CASE (uniform)
538 rng_stream%distribution_type = uniform
539 CASE DEFAULT
540 cpabort("Invalid distribution type specified")
541 END SELECT
542 ELSE IF (PRESENT(last_rng_stream)) THEN
543 rng_stream%distribution_type = last_rng_stream%distribution_type
544 END IF
545
546 IF (PRESENT(antithetic)) THEN
547 rng_stream%antithetic = antithetic
548 ELSE IF (PRESENT(last_rng_stream)) THEN
549 rng_stream%antithetic = last_rng_stream%antithetic
550 END IF
551
552 IF (PRESENT(extended_precision)) THEN
553 rng_stream%extended_precision = extended_precision
554 ELSE IF (PRESENT(last_rng_stream)) THEN
555 rng_stream%extended_precision = last_rng_stream%extended_precision
556 END IF
557 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 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 rng_stream%name, &
570 rng_stream%distribution_type, &
571 rng_stream%antithetic, &
572 rng_stream%extended_precision, &
573 rng_stream%buffer_filled, &
574 rng_stream%buffer, &
575 rng_stream%cg, &
576 rng_stream%bg, &
577 rng_stream%ig
578 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 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 rng_record = " "
590 WRITE (unit=rng_record, fmt=rng_record_format) &
591 self%name, &
592 self%distribution_type, &
593 self%antithetic, &
594 self%extended_precision, &
595 self%buffer_filled, &
596 self%buffer, &
597 self%cg, &
598 self%bg, &
599 self%ig
600 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 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 IF (PRESENT(name)) name = self%name
634 IF (PRESENT(distribution_type)) THEN
635 distribution_type = self%distribution_type
636 END IF
637 IF (PRESENT(bg)) bg = self%bg
638 IF (PRESENT(cg)) cg = self%cg
639 IF (PRESENT(ig)) ig = self%ig
640 IF (PRESENT(antithetic)) antithetic = self%antithetic
641 IF (PRESENT(extended_precision)) THEN
642 extended_precision = self%extended_precision
643 END IF
644 IF (PRESENT(buffer)) buffer = self%buffer
645 IF (PRESENT(buffer_filled)) buffer_filled = self%buffer_filled
646 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 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 DO i = 1, 3
663 CALL mat_vec_mod_m(a, b(:, i), c(:, i), m)
664 END DO
665
666 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 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 w = a
687
688 b(1, 1) = 1.0_dp
689 b(2, 1) = 0.0_dp
690 b(3, 1) = 0.0_dp
691 b(1, 2) = 0.0_dp
692 b(2, 2) = 1.0_dp
693 b(3, 2) = 0.0_dp
694 b(1, 3) = 0.0_dp
695 b(2, 3) = 0.0_dp
696 b(3, 3) = 1.0_dp
697
698 ! Compute b = MODULO(a**n,m) using the binary decomposition of n
699
700 i = n
701
702 DO
703 IF (modulo(i, 2) /= 0) THEN
704 u = w
705 v = b
706 CALL mat_mat_mod_m(u, v, b, m)
707 END IF
708 i = i/2
709 IF (i == 0) EXIT
710 u = w
711 v = w
712 CALL mat_mat_mod_m(u, v, w, m)
713 END DO
714 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 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 b = a
733
734 DO i = 1, e
735 u = b
736 v = b
737 CALL mat_mat_mod_m(u, v, b, m)
738 END DO
739
740 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 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 v = 0.0_dp
759
760 DO i = 1, 3
761 DO j = 1, 3
762 a2 = a(i, j)
763 c = v(i)
764 v(i) = a2*s(j) + c
765 IF ((v(i) >= two53) .OR. (v(i) <= -two53)) THEN
766 a1 = int(a2/two17)
767 a2 = a2 - a1*two17
768 v(i) = a1*s(j)
769 a1 = int(v(i)/m)
770 v(i) = v(i) - a1*m
771 v(i) = v(i)*two17 + a2*s(j) + c
772 END IF
773 a1 = int(v(i)/m)
774 v(i) = v(i) - a1*m
775 IF (v(i) < 0.0_dp) v(i) = v(i) + m
776 END DO
777 END DO
778
779 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 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 cpassert(self%distribution_type == uniform)
796
797 r = self%next_real()
798 u = low + int(r*real(high - low + 1, dp))
799 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 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 SELECT CASE (self%distribution_type)
815 CASE (gaussian)
816 var = 1.0_dp
817 IF (PRESENT(variance)) var = variance
818 ! take the random number from the buffer, if the buffer is filled
819 IF (self%buffer_filled) THEN
820 u = sqrt(var)*self%buffer
821 self%buffer_filled = .false.
822 ELSE
823 DO
824 IF (self%extended_precision) THEN
825 u1 = 2.0_dp*rn53(self) - 1.0_dp
826 u2 = 2.0_dp*rn53(self) - 1.0_dp
827 ELSE
828 u1 = 2.0_dp*rn32(self) - 1.0_dp
829 u2 = 2.0_dp*rn32(self) - 1.0_dp
830 END IF
831 r = u1*u1 + u2*u2
832 IF ((r > 0.0_dp) .AND. (r < 1.0_dp)) EXIT
833 END DO
834 ! Box-Muller transformation
835 f = sqrt(-2.0_dp*log(r)/r)
836 u = sqrt(var)*f*u1
837 ! save the second random number for the next call
838 self%buffer = f*u2
839 self%buffer_filled = .true.
840 END IF
841 CASE (uniform)
842 IF (self%extended_precision) THEN
843 u = rn53(self)
844 ELSE
845 u = rn32(self)
846 END IF
847 END SELECT
848 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 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 IF (PRESENT(seed)) THEN
861 CALL check_seed(seed)
862 CALL mat_vec_mod_m(a1p127, seed(:, 1), next_seed(:, 1), m1)
863 CALL mat_vec_mod_m(a2p127, seed(:, 2), next_seed(:, 2), m2)
864 ELSE
865 next_seed = 12345.0_dp ! default seed
866 END IF
867
868 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 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 DO i = 1, SIZE(array)
882 array(i) = self%next()
883 END DO
884 END SUBROUTINE fill_1
885
886! **************************************************************************************************
887!> \brief ...
888!> \param self ...
889!> \param array ...
890! **************************************************************************************************
891 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 DO j = 1, SIZE(array, 2)
898 DO i = 1, SIZE(array, 1)
899 array(i, j) = self%next()
900 END DO
901 END DO
902 END SUBROUTINE fill_2
903
904! **************************************************************************************************
905!> \brief ...
906!> \param self ...
907!> \param array ...
908! **************************************************************************************************
909 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 DO k = 1, SIZE(array, 3)
916 DO j = 1, SIZE(array, 2)
917 DO i = 1, SIZE(array, 1)
918 array(i, j, k) = self%next()
919 END DO
920 END DO
921 END DO
922 END SUBROUTINE fill_3
923
924! **************************************************************************************************
925!> \brief Reset a random number stream to its initial state.
926!> \param self ...
927! **************************************************************************************************
928 SUBROUTINE reset(self)
929 CLASS(rng_stream_type), INTENT(INOUT) :: self
930
931 self%cg = self%ig
932 self%bg = self%ig
933 END SUBROUTINE reset
934
935! **************************************************************************************************
936!> \brief Reset a random number stream to the beginning of its current substream.
937!> \param self ...
938! **************************************************************************************************
939 SUBROUTINE reset_to_substream(self)
940 CLASS(rng_stream_type), INTENT(INOUT) :: self
941
942 self%cg = self%bg
943 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 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 u = 0.0_dp
955
956 CALL mat_vec_mod_m(a1p76, self%bg(:, 1), u(:, 1), m1)
957 CALL mat_vec_mod_m(a2p76, self%bg(:, 2), u(:, 2), m2)
958
959 self%bg = u
960 self%cg = u
961 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 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 p1 = a12*rng_stream%cg(2, 1) - a13n*rng_stream%cg(1, 1)
978 k = int(p1/m1)
979 p1 = p1 - k*m1
980 IF (p1 < 0.0_dp) p1 = p1 + m1
981 rng_stream%cg(1, 1) = rng_stream%cg(2, 1)
982 rng_stream%cg(2, 1) = rng_stream%cg(3, 1)
983 rng_stream%cg(3, 1) = p1
984
985 ! Component 2
986
987 p2 = a21*rng_stream%cg(3, 2) - a23n*rng_stream%cg(1, 2)
988 k = int(p2/m2)
989 p2 = p2 - k*m2
990 IF (p2 < 0.0_dp) p2 = p2 + m2
991 rng_stream%cg(1, 2) = rng_stream%cg(2, 2)
992 rng_stream%cg(2, 2) = rng_stream%cg(3, 2)
993 rng_stream%cg(3, 2) = p2
994
995 ! Combination
996
997 IF (p1 > p2) THEN
998 u = (p1 - p2)*norm
999 ELSE
1000 u = (p1 - p2 + m1)*norm
1001 END IF
1002
1003 IF (rng_stream%antithetic) u = 1.0_dp - u
1004
1005 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 FUNCTION rn53(rng_stream) RESULT(u)
1013 TYPE(rng_stream_type) :: rng_stream
1014 REAL(kind=dp) :: u
1015
1016 u = rn32(rng_stream)
1017
1018 ! Note: rn32 returns 1 - u in the antithetic case
1019
1020 IF (rng_stream%antithetic) THEN
1021 u = u + (rn32(rng_stream) - 1.0_dp)*fact
1022 IF (u < 0.0_dp) u = u + 1.0_dp
1023 ELSE
1024 u = u + rn32(rng_stream)*fact
1025 IF (u >= 1.0_dp) u = u - 1.0_dp
1026 END IF
1027 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 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 IF (PRESENT(name)) self%name = name
1062 IF (PRESENT(distribution_type)) THEN
1063 self%distribution_type = distribution_type
1064 END IF
1065 IF (PRESENT(bg)) self%bg = bg
1066 IF (PRESENT(cg)) self%cg = cg
1067 IF (PRESENT(ig)) self%ig = ig
1068 IF (PRESENT(seed)) THEN
1069 ! Sets the initial seed of the stream to seed
1070 ! NOTE: The use of this method is discouraged
1071 CALL check_seed(seed)
1072 self%ig = seed
1073 self%cg = seed
1074 self%bg = seed
1075 END IF
1076 IF (PRESENT(antithetic)) self%antithetic = antithetic
1077 IF (PRESENT(extended_precision)) THEN
1078 self%extended_precision = extended_precision
1079 END IF
1080 IF (PRESENT(buffer)) self%buffer = buffer
1081 IF (PRESENT(buffer_filled)) self%buffer_filled = buffer_filled
1082 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 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 "GENERATOR"
1099
1100 fmtstr = "(/,T4,A,/,/,(2X,3F14.1))"
1101
1102 WRITE (unit=output_unit, fmt=fmtstr) &
1103 "A1", ((a1p0(i, j), j=1, 3), i=1, 3)
1104
1105 WRITE (unit=output_unit, fmt=fmtstr) &
1106 "A2", ((a2p0(i, j), j=1, 3), i=1, 3)
1107
1108 WRITE (unit=output_unit, fmt=fmtstr) &
1109 "A1**(2**76)", ((a1p76(i, j), j=1, 3), i=1, 3)
1110
1111 WRITE (unit=output_unit, fmt=fmtstr) &
1112 "A2**(2**76)", ((a2p76(i, j), j=1, 3), i=1, 3)
1113
1114 WRITE (unit=output_unit, fmt=fmtstr) &
1115 "A1**(2**127)", ((a1p127(i, j), j=1, 3), i=1, 3)
1116
1117 WRITE (unit=output_unit, fmt=fmtstr) &
1118 "A2**(2**127)", ((a2p127(i, j), j=1, 3), i=1, 3)
1119
1120 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 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 my_write_all = .false.
1136
1137 IF (PRESENT(write_all)) THEN
1138 my_write_all = write_all
1139 END IF
1140
1141 WRITE (unit=output_unit, fmt="(/,T2,A,/)") &
1142 "Random number stream <"//trim(self%name)//">:"
1143
1144 SELECT CASE (self%distribution_type)
1145 CASE (gaussian)
1146 WRITE (unit=output_unit, fmt="(T4,A)") &
1147 "Distribution type: "// &
1148 "Normal Gaussian distribution with zero mean"
1149 CASE (uniform)
1150 WRITE (unit=output_unit, fmt="(T4,A)") &
1151 "Distribution type: "// &
1152 "Uniform distribution [0,1] with 1/2 mean"
1153 END SELECT
1154
1155 IF (self%antithetic) THEN
1156 WRITE (unit=output_unit, fmt="(T4,A)") "Antithetic: yes"
1157 ELSE
1158 WRITE (unit=output_unit, fmt="(T4,A)") "Antithetic: no"
1159 END IF
1160
1161 IF (self%extended_precision) THEN
1162 WRITE (unit=output_unit, fmt="(T4,A)") "Precision: 53 Bit"
1163 ELSE
1164 WRITE (unit=output_unit, fmt="(T4,A)") "Precision: 32 Bit"
1165 END IF
1166
1167 IF (my_write_all) THEN
1168
1169 WRITE (unit=output_unit, fmt="(/,T4,A,/,/,(T4,A,3F20.1))") &
1170 "Initial state of the stream:", &
1171 "Component 1:", self%ig(:, 1), &
1172 "Component 2:", self%ig(:, 2)
1173
1174 WRITE (unit=output_unit, fmt="(/,T4,A,/,/,(T4,A,3F20.1))") &
1175 "Initial state of the current substream:", &
1176 "Component 1:", self%bg(:, 1), &
1177 "Component 2:", self%bg(:, 2)
1178
1179 END IF
1180
1181 WRITE (unit=output_unit, fmt="(/,T4,A,/,/,(T4,A,3F20.1))") &
1182 "Current state of the stream:", &
1183 "Component 1:", self%cg(:, 1), &
1184 "Component 2:", self%cg(:, 2)
1185 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 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 DO idxa = ubound(arr, 1), lbound(arr, 1) + 1, -1
1199 idxb = self%next(lbound(arr, 1), idxa)
1200 tmp = arr(idxa)
1201 arr(idxa) = arr(idxb)
1202 arr(idxb) = tmp
1203 END DO
1204 END SUBROUTINE shuffle
1205
1206END 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:125
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:565
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:856
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:1089
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:3191
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