(git:0de0cc2)
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 ! **************************************************************************************************
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, &
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  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 
1202 END MODULE parallel_rng_types
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
Defines the basic variable types.
Definition: kinds.F:23
integer, parameter, public dp
Definition: kinds.F:34
integer, parameter, public default_string_length
Definition: kinds.F:57
Parallel (pseudo)random number generator (RNG) for multiple streams and substreams of random numbers.
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).
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.
integer, parameter, public rng_name_length
integer, parameter, public rng_record_length
integer, parameter, public uniform
subroutine, public check_rng(output_unit, ionode)
...
subroutine, public write_rng_matrices(output_unit)
Write the transformation matrices of the two MRG components (raised to the specified output)
integer, parameter, public gaussian
Utilities for string manipulations.
subroutine, public compress(string, full)
Eliminate multiple space characters in a string. If full is .TRUE., then all spaces are eliminated.