56 #include "../base/base_uses.f90"
64 CHARACTER(LEN=*),
PARAMETER,
PRIVATE :: rng_record_format =
"(A40,I2,3L2,ES25.16,18F20.1)"
76 REAL(kind=
dp),
PARAMETER :: norm = 2.328306549295727688e-10_dp, &
77 m1 = 4294967087.0_dp, &
78 m2 = 4294944443.0_dp, &
82 a23n = 1370589.0_dp, &
83 two17 = 131072.0_dp, &
84 two53 = 9007199254740992.0_dp, &
85 fact = 5.9604644775390625e-8_dp
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 &
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 &
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 &
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 &
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 &
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 &
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 &
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 &
161 CHARACTER(LEN=rng_name_length) :: name =
""
162 INTEGER :: distribution_type =
uniform
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.
167 REAL(kind=
dp) :: buffer = 0.0_dp
168 LOGICAL :: buffer_filled = .false.
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
176 PROCEDURE, pass(self) :: next_int
177 PROCEDURE, pass(self) :: next_real
178 generic,
PUBLIC :: next => next_int, next_real
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
191 INTERFACE rng_stream_type
192 MODULE PROCEDURE :: rng_stream_constructor
195 TYPE rng_stream_p_type
196 TYPE(rng_stream_type),
POINTER :: stream => null()
197 END TYPE rng_stream_p_type
208 PUBLIC :: rng_stream_p_type, &
227 SUBROUTINE advance(self, e, c)
228 CLASS(rng_stream_type),
INTENT(INOUT) :: self
229 INTEGER,
INTENT(IN) :: e, c
231 REAL(kind=
dp),
DIMENSION(3, 2) :: x
232 REAL(kind=
dp),
DIMENSION(3, 3) :: u1, u2, v1, v2, w1, w2
242 CALL mat_two_pow_mod_m(a1p0, u1, m1, e)
243 CALL mat_two_pow_mod_m(a2p0, u2, m2, e)
245 CALL mat_two_pow_mod_m(inv_a1, u1, m1, -e)
246 CALL mat_two_pow_mod_m(inv_a2, u2, m2, -e)
250 CALL mat_pow_mod_m(a1p0, v1, m1, c)
251 CALL mat_pow_mod_m(a2p0, v2, m2, c)
253 CALL mat_pow_mod_m(inv_a1, v1, m1, -c)
254 CALL mat_pow_mod_m(inv_a2, v2, m2, -c)
261 CALL mat_mat_mod_m(u1, v1, w1, m1)
262 CALL mat_mat_mod_m(u2, v2, w2, m2)
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)
271 END SUBROUTINE advance
282 INTEGER,
INTENT(IN) :: output_unit
283 LOGICAL,
INTENT(IN) :: ionode
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, &
296 g1 = rng_stream_type(
"g1")
297 g2 = rng_stream_type(
"g2", g1)
298 g3 = rng_stream_type(
"g3", g2)
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)
309 sum = g2%next() + g3%next()
311 CALL g1%advance(5, 3)
312 sum = sum + g1%next()
316 CALL g1%advance(0, 1)
318 sum = sum + g1%next()
324 sumi = sumi + g1%next(1, 10)
326 sum = sum + sumi/100.0_dp
330 sum3 = sum3 + g3%next()
332 sum = sum + sum3/10.0_dp
336 sum = sum + g3%next()
341 CALL g3%reset_to_next_substream()
344 sum = sum + g3%next()
347 CALL g3%reset_to_substream()
349 sum = sum + g3%next()
352 CALL g2%reset_to_next_substream()
355 sum3 = sum3 + g2%next()
357 sum = sum + sum3/10000.0_dp
359 CALL g3%set(antithetic=.true.)
362 sum3 = sum3 + g3%next()
364 sum = sum + sum3/10000.0_dp
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
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)
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)
397 sum = sum + poisson%next() + laplace%next() + galois%next() + cantor%next()
399 CALL galois%advance(-127, 0)
400 sum = sum + galois%next()
402 CALL galois%reset_to_next_substream()
403 CALL galois%set(extended_precision=.true.)
406 sum3 = sum3 + galois%next()
408 sum = sum + sum3/10000.0_dp
410 CALL galois%set(antithetic=.true.)
413 sum3 = sum3 + galois%next()
415 sum = sum + sum3/10000.0_dp
416 CALL galois%set(antithetic=.false.)
418 CALL galois%set(extended_precision=.false.)
419 sum = sum + poisson%next() + laplace%next() + galois%next() + cantor%next()
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
440 SUBROUTINE check_seed(seed)
441 REAL(kind=
dp),
DIMENSION(3, 2),
INTENT(IN) :: seed
443 CHARACTER(LEN=*),
PARAMETER :: fmtstr =
"(A,I1,A,ES23.14,A,ES23.14)"
445 CHARACTER(LEN=default_string_length) :: message
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
458 IF (seed(i, 1) >= m1)
THEN
459 WRITE (unit=message, fmt=fmtstr) &
460 "seed(", i,
",1) = ", seed(i, 1),
" >= ", m1
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
473 IF (seed(i, 2) >= m2)
THEN
474 WRITE (unit=message, fmt=fmtstr) &
475 "seed(", i,
",2) = ", seed(i, 2),
" >= ", m2
484 IF (all(seed(:, 1) < 1.0_dp))
THEN
485 cpabort(
"First seed = 0")
488 IF (all(seed(:, 2) < 1.0_dp))
THEN
489 cpabort(
"Second seed = 0")
492 END SUBROUTINE check_seed
504 FUNCTION rng_stream_constructor(name, last_rng_stream, distribution_type, seed, antithetic, extended_precision) &
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), &
512 LOGICAL,
INTENT(IN),
OPTIONAL :: antithetic, extended_precision
513 TYPE(rng_stream_type) :: rng_stream
516 cpabort(
"given random number generator name is too long")
518 rng_stream%name = trim(name)
520 IF (
PRESENT(seed))
THEN
521 CALL check_seed(seed)
523 ELSE IF (
PRESENT(last_rng_stream))
THEN
529 rng_stream%cg = rng_stream%ig
530 rng_stream%bg = rng_stream%ig
532 IF (
PRESENT(distribution_type))
THEN
533 SELECT CASE (distribution_type)
535 rng_stream%distribution_type =
gaussian
537 rng_stream%distribution_type =
uniform
539 cpabort(
"Invalid distribution type specified")
541 ELSE IF (
PRESENT(last_rng_stream))
THEN
542 rng_stream%distribution_type = last_rng_stream%distribution_type
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
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
556 END FUNCTION rng_stream_constructor
564 CHARACTER(LEN=rng_record_length),
INTENT(IN) :: rng_record
565 TYPE(rng_stream_type) :: rng_stream
567 READ (unit=rng_record, fmt=rng_record_format) &
569 rng_stream%distribution_type, &
570 rng_stream%antithetic, &
571 rng_stream%extended_precision, &
572 rng_stream%buffer_filled, &
584 SUBROUTINE dump(self, rng_record)
585 CLASS(rng_stream_type),
INTENT(IN) :: self
586 CHARACTER(LEN=rng_record_length),
INTENT(OUT) :: rng_record
589 WRITE (unit=rng_record, fmt=rng_record_format) &
591 self%distribution_type, &
593 self%extended_precision, &
594 self%buffer_filled, &
619 SUBROUTINE get(self, name, distribution_type, bg, cg, ig, &
620 antithetic, extended_precision, &
621 buffer, buffer_filled)
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
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
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
660 CALL mat_vec_mod_m(a, b(:, i), c(:, i), m)
663 END SUBROUTINE mat_mat_mod_m
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
679 REAL(kind=
dp),
DIMENSION(3, 3) :: u, v, w
700 IF (
modulo(i, 2) /= 0)
THEN
703 CALL mat_mat_mod_m(u, v, b, m)
709 CALL mat_mat_mod_m(u, v, w, m)
711 END SUBROUTINE mat_pow_mod_m
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
727 REAL(kind=
dp),
DIMENSION(3, 3) :: u, v
734 CALL mat_mat_mod_m(u, v, b, m)
737 END SUBROUTINE mat_two_pow_mod_m
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
753 REAL(kind=
dp) :: a1, a2, c
762 IF ((v(i) >= two53) .OR. (v(i) <= -two53))
THEN
768 v(i) = v(i)*two17 + a2*s(j) + c
772 IF (v(i) < 0.0_dp) v(i) = v(i) + m
776 END SUBROUTINE mat_vec_mod_m
785 FUNCTION next_int(self, low, high)
RESULT(u)
786 CLASS(rng_stream_type),
INTENT(INOUT) :: self
787 INTEGER,
INTENT(IN) :: low, high
792 cpassert(self%distribution_type ==
uniform)
795 u = low + int(r*real(high - low + 1,
dp))
796 END FUNCTION next_int
804 FUNCTION next_real(self, variance)
RESULT(u)
805 CLASS(rng_stream_type),
INTENT(INOUT) :: self
806 REAL(kind=
dp),
INTENT(IN),
OPTIONAL :: variance
809 REAL(kind=
dp) :: f, r, u1, u2, var
811 SELECT CASE (self%distribution_type)
814 IF (
PRESENT(variance)) var = variance
816 IF (self%buffer_filled)
THEN
817 u = sqrt(var)*self%buffer
818 self%buffer_filled = .false.
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
825 u1 = 2.0_dp*rn32(self) - 1.0_dp
826 u2 = 2.0_dp*rn32(self) - 1.0_dp
829 IF ((r > 0.0_dp) .AND. (r < 1.0_dp))
EXIT
832 f = sqrt(-2.0_dp*log(r)/r)
836 self%buffer_filled = .true.
839 IF (self%extended_precision)
THEN
845 END FUNCTION next_real
853 REAL(kind=
dp),
DIMENSION(3, 2),
INTENT(IN), &
855 REAL(kind=
dp),
DIMENSION(3, 2) :: next_seed
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)
862 next_seed = 12345.0_dp
872 SUBROUTINE fill_1(self, array)
873 CLASS(rng_stream_type),
INTENT(INOUT) :: self
874 REAL(kind=
dp),
DIMENSION(:),
INTENT(OUT) :: array
878 DO i = 1,
SIZE(array)
879 array(i) = self%next()
881 END SUBROUTINE fill_1
888 SUBROUTINE fill_2(self, array)
889 CLASS(rng_stream_type),
INTENT(INOUT) :: self
890 REAL(kind=
dp),
DIMENSION(:, :),
INTENT(OUT) :: array
894 DO j = 1,
SIZE(array, 2)
895 DO i = 1,
SIZE(array, 1)
896 array(i, j) = self%next()
899 END SUBROUTINE fill_2
906 SUBROUTINE fill_3(self, array)
907 CLASS(rng_stream_type),
INTENT(INOUT) :: self
908 REAL(kind=
dp),
DIMENSION(:, :, :),
INTENT(OUT) :: array
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()
919 END SUBROUTINE fill_3
925 SUBROUTINE reset(self)
926 CLASS(rng_stream_type),
INTENT(INOUT) :: self
936 SUBROUTINE reset_to_substream(self)
937 CLASS(rng_stream_type),
INTENT(INOUT) :: self
940 END SUBROUTINE reset_to_substream
946 SUBROUTINE reset_to_next_substream(self)
947 CLASS(rng_stream_type),
INTENT(INOUT) :: self
949 REAL(kind=
dp),
DIMENSION(3, 2) :: u
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)
958 END SUBROUTINE reset_to_next_substream
965 FUNCTION rn32(rng_stream)
RESULT(u)
966 TYPE(rng_stream_type) :: rng_stream
970 REAL(kind=
dp) :: p1, p2
974 p1 = a12*rng_stream%cg(2, 1) - a13n*rng_stream%cg(1, 1)
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
984 p2 = a21*rng_stream%cg(3, 2) - a23n*rng_stream%cg(1, 2)
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
997 u = (p1 - p2 + m1)*norm
1000 IF (rng_stream%antithetic) u = 1.0_dp - u
1009 FUNCTION rn53(rng_stream)
RESULT(u)
1010 TYPE(rng_stream_type) :: rng_stream
1013 u = rn32(rng_stream)
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
1021 u = u + rn32(rng_stream)*fact
1022 IF (u >= 1.0_dp) u = u - 1.0_dp
1043 SUBROUTINE set(self, name, distribution_type, bg, cg, ig, &
1044 seed, antithetic, extended_precision, &
1045 buffer, buffer_filled)
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
1058 IF (
PRESENT(name)) self%name = name
1059 IF (
PRESENT(distribution_type))
THEN
1060 self%distribution_type = distribution_type
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
1068 CALL check_seed(seed)
1073 IF (
PRESENT(antithetic)) self%antithetic = antithetic
1074 IF (
PRESENT(extended_precision))
THEN
1075 self%extended_precision = extended_precision
1077 IF (
PRESENT(buffer)) self%buffer = buffer
1078 IF (
PRESENT(buffer_filled)) self%buffer_filled = buffer_filled
1086 INTEGER,
INTENT(IN) :: output_unit
1088 CHARACTER(LEN=40) :: fmtstr
1093 WRITE (unit=output_unit, fmt=
"(/,T2,A)") &
1094 "TRANSFORMATION MATRICES FOR THE PARALLEL (PSEUDO)RANDOM NUMBER "// &
1097 fmtstr =
"(/,T4,A,/,/,(2X,3F14.1))"
1099 WRITE (unit=output_unit, fmt=fmtstr) &
1100 "A1", ((a1p0(i, j), j=1, 3), i=1, 3)
1102 WRITE (unit=output_unit, fmt=fmtstr) &
1103 "A2", ((a2p0(i, j), j=1, 3), i=1, 3)
1105 WRITE (unit=output_unit, fmt=fmtstr) &
1106 "A1**(2**76)", ((a1p76(i, j), j=1, 3), i=1, 3)
1108 WRITE (unit=output_unit, fmt=fmtstr) &
1109 "A2**(2**76)", ((a2p76(i, j), j=1, 3), i=1, 3)
1111 WRITE (unit=output_unit, fmt=fmtstr) &
1112 "A1**(2**127)", ((a1p127(i, j), j=1, 3), i=1, 3)
1114 WRITE (unit=output_unit, fmt=fmtstr) &
1115 "A2**(2**127)", ((a2p127(i, j), j=1, 3), i=1, 3)
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
1130 LOGICAL :: my_write_all
1132 my_write_all = .false.
1134 IF (
PRESENT(write_all)) &
1135 my_write_all = write_all
1137 WRITE (unit=output_unit, fmt=
"(/,T2,A,/)") &
1138 "Random number stream <"//trim(self%name)//
">:"
1140 SELECT CASE (self%distribution_type)
1142 WRITE (unit=output_unit, fmt=
"(T4,A)") &
1143 "Distribution type: "// &
1144 "Normal Gaussian distribution with zero mean"
1146 WRITE (unit=output_unit, fmt=
"(T4,A)") &
1147 "Distribution type: "// &
1148 "Uniform distribution [0,1] with 1/2 mean"
1151 IF (self%antithetic)
THEN
1152 WRITE (unit=output_unit, fmt=
"(T4,A)")
"Antithetic: yes"
1154 WRITE (unit=output_unit, fmt=
"(T4,A)")
"Antithetic: no"
1157 IF (self%extended_precision)
THEN
1158 WRITE (unit=output_unit, fmt=
"(T4,A)")
"Precision: 53 Bit"
1160 WRITE (unit=output_unit, fmt=
"(T4,A)")
"Precision: 32 Bit"
1163 IF (my_write_all)
THEN
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)
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)
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
1188 SUBROUTINE shuffle(self, arr)
1189 CLASS(rng_stream_type),
INTENT(INOUT) :: self
1190 INTEGER,
DIMENSION(:),
INTENT(INOUT) :: arr
1192 INTEGER :: idxa, idxb, tmp
1194 DO idxa = ubound(arr, 1), lbound(arr, 1) + 1, -1
1195 idxb = self%next(lbound(arr, 1), idxa)
1197 arr(idxa) = arr(idxb)
static GRID_HOST_DEVICE int modulo(int a, int m)
Equivalent of Fortran's MODULO, which always return a positive number. https://gcc....
Defines the basic variable types.
integer, parameter, public dp
integer, parameter, public default_string_length
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.