|
111
|
1 ! { dg-do compile }
|
|
|
2 ! { dg-options "-O2 -fdump-tree-optimized" }
|
|
|
3 !
|
|
|
4 ! PR fortran/54556
|
|
|
5 !
|
|
|
6 ! Contributed by Joost VandeVondele
|
|
|
7 !
|
|
|
8 MODULE parallel_rng_types
|
|
|
9
|
|
|
10 IMPLICIT NONE
|
|
|
11
|
|
|
12 ! Global parameters in this module
|
|
|
13 INTEGER, PARAMETER :: dp=8
|
|
|
14
|
|
|
15 TYPE rng_stream_type
|
|
|
16 PRIVATE
|
|
|
17 CHARACTER(LEN=40) :: name
|
|
|
18 INTEGER :: distribution_type
|
|
|
19 REAL(KIND=dp), DIMENSION(3,2) :: bg,cg,ig
|
|
|
20 LOGICAL :: antithetic,extended_precision
|
|
|
21 REAL(KIND=dp) :: buffer
|
|
|
22 LOGICAL :: buffer_filled
|
|
|
23 END TYPE rng_stream_type
|
|
|
24
|
|
|
25 REAL(KIND=dp), DIMENSION(3,3) :: a1p0,a1p76,a1p127,&
|
|
|
26 a2p0,a2p76,a2p127,&
|
|
|
27 inv_a1,inv_a2
|
|
|
28
|
|
|
29 INTEGER, PARAMETER :: GAUSSIAN = 1,&
|
|
|
30 UNIFORM = 2
|
|
|
31
|
|
|
32 REAL(KIND=dp), PARAMETER :: norm = 2.328306549295727688e-10_dp,&
|
|
|
33 m1 = 4294967087.0_dp,&
|
|
|
34 m2 = 4294944443.0_dp,&
|
|
|
35 a12 = 1403580.0_dp,&
|
|
|
36 a13n = 810728.0_dp,&
|
|
|
37 a21 = 527612.0_dp,&
|
|
|
38 a23n = 1370589.0_dp,&
|
|
|
39 two17 = 131072.0_dp,& ! 2**17
|
|
|
40 two53 = 9007199254740992.0_dp,& ! 2**53
|
|
|
41 fact = 5.9604644775390625e-8_dp ! 1/2**24
|
|
|
42
|
|
|
43
|
|
|
44 CONTAINS
|
|
|
45
|
|
|
46 FUNCTION rn32(rng_stream) RESULT(u)
|
|
|
47
|
|
|
48 TYPE(rng_stream_type), POINTER :: rng_stream
|
|
|
49 REAL(KIND=dp) :: u
|
|
|
50
|
|
|
51 INTEGER :: k
|
|
|
52 REAL(KIND=dp) :: p1, p2
|
|
|
53
|
|
|
54 ! -------------------------------------------------------------------------
|
|
|
55 ! Component 1
|
|
|
56
|
|
|
57 p1 = a12*rng_stream%cg(2,1) - a13n*rng_stream%cg(1,1)
|
|
|
58 k = INT(p1/m1)
|
|
|
59 p1 = p1 - k*m1
|
|
|
60 IF (p1 < 0.0_dp) p1 = p1 + m1
|
|
|
61 rng_stream%cg(1,1) = rng_stream%cg(2,1)
|
|
|
62 rng_stream%cg(2,1) = rng_stream%cg(3,1)
|
|
|
63 rng_stream%cg(3,1) = p1
|
|
|
64
|
|
|
65 ! Component 2
|
|
|
66
|
|
|
67 p2 = a21*rng_stream%cg(3,2) - a23n*rng_stream%cg(1,2)
|
|
|
68 k = INT(p2/m2)
|
|
|
69 p2 = p2 - k*m2
|
|
|
70 IF (p2 < 0.0_dp) p2 = p2 + m2
|
|
|
71 rng_stream%cg(1,2) = rng_stream%cg(2,2)
|
|
|
72 rng_stream%cg(2,2) = rng_stream%cg(3,2)
|
|
|
73 rng_stream%cg(3,2) = p2
|
|
|
74
|
|
|
75 ! Combination
|
|
|
76
|
|
|
77 IF (p1 > p2) THEN
|
|
|
78 u = (p1 - p2)*norm
|
|
|
79 ELSE
|
|
|
80 u = (p1 - p2 + m1)*norm
|
|
|
81 END IF
|
|
|
82
|
|
|
83 IF (rng_stream%antithetic) u = 1.0_dp - u
|
|
|
84
|
|
|
85 END FUNCTION rn32
|
|
|
86
|
|
|
87 ! *****************************************************************************
|
|
|
88 FUNCTION rn53(rng_stream) RESULT(u)
|
|
|
89
|
|
|
90 TYPE(rng_stream_type), POINTER :: rng_stream
|
|
|
91 REAL(KIND=dp) :: u
|
|
|
92
|
|
|
93 u = rn32(rng_stream)
|
|
|
94
|
|
|
95 IF (rng_stream%antithetic) THEN
|
|
|
96 u = u + (rn32(rng_stream) - 1.0_dp)*fact
|
|
|
97 IF (u < 0.0_dp) u = u + 1.0_dp
|
|
|
98 ELSE
|
|
|
99 u = u + rn32(rng_stream)*fact
|
|
|
100 IF (u >= 1.0_dp) u = u - 1.0_dp
|
|
|
101 END IF
|
|
|
102
|
|
|
103 END FUNCTION rn53
|
|
|
104
|
|
|
105 END MODULE
|
|
|
106
|
|
|
107 ! { dg-final { scan-module-absence "parallel_rng_types" "IMPLICIT_PURE" } }
|
|
|
108 ! { dg-final { scan-tree-dump-times "rn32 \\(rng_stream" 3 "optimized" } }
|
