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