1! { dg-do compile }
2! { dg-options "-O2 -fdump-tree-optimized" }
3!
4! PR fortran/54556
5!
6! Contributed by Joost VandeVondele
7!
8MODULE 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
44CONTAINS
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
105END 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" } }
109