Line data Source code
1 : ! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2 : ! Copyright (c) 2015, Regents of the University of Colorado
3 : ! All rights reserved.
4 : !
5 : ! Redistribution and use in source and binary forms, with or without modification, are
6 : ! permitted provided that the following conditions are met:
7 : !
8 : ! 1. Redistributions of source code must retain the above copyright notice, this list of
9 : ! conditions and the following disclaimer.
10 : !
11 : ! 2. Redistributions in binary form must reproduce the above copyright notice, this list
12 : ! of conditions and the following disclaimer in the documentation and/or other
13 : ! materials provided with the distribution.
14 : !
15 : ! 3. Neither the name of the copyright holder nor the names of its contributors may be
16 : ! used to endorse or promote products derived from this software without specific prior
17 : ! written permission.
18 : !
19 : ! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY
20 : ! EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
21 : ! MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL
22 : ! THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
23 : ! SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT
24 : ! OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
25 : ! INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
26 : ! LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
27 : ! OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
28 : !
29 : ! History:
30 : ! May 2015- D. Swales - Original version
31 : !
32 : ! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
33 : MODULE mod_rng
34 : USE cosp_kinds, ONLY: dp, sp, wp
35 : IMPLICIT NONE
36 :
37 : INTEGER, parameter :: ki9 = selected_int_kind(R=9)
38 : integer :: testInt
39 :
40 : TYPE rng_state
41 : INTEGER(ki9) :: seed ! 32 bit integer
42 : END TYPE rng_state
43 :
44 : INTERFACE init_rng
45 : MODULE PROCEDURE init_rng_1, init_rng_n
46 : END INTERFACE init_rng
47 :
48 : INTERFACE get_rng
49 : MODULE PROCEDURE get_rng_1, get_rng_n, get_rng_v
50 : END INTERFACE get_rng
51 :
52 : CONTAINS
53 :
54 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
55 : ! Set single seed
56 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
57 0 : SUBROUTINE init_rng_1(s, seed_in)
58 : TYPE(rng_state), INTENT(INOUT) :: s
59 : INTEGER, INTENT(IN ) :: seed_in
60 0 : s%seed = seed_in
61 0 : END SUBROUTINE init_rng_1
62 :
63 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
64 : ! Set vector of seeds
65 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
66 9288 : SUBROUTINE init_rng_n(s, seed_in)
67 : TYPE(rng_state), DIMENSION(:), INTENT(INOUT) :: s
68 : INTEGER, DIMENSION(:), INTENT(IN ) :: seed_in
69 :
70 : INTEGER :: i
71 155088 : DO i = 1, SIZE(seed_in)
72 155088 : s(i)%seed = seed_in(i)
73 : END DO
74 9288 : END SUBROUTINE init_rng_n
75 :
76 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
77 : ! Create single random number from seed
78 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
79 123930000 : FUNCTION get_rng_1(s)
80 : TYPE(rng_state), INTENT(INOUT) :: s
81 : REAL(WP) :: get_rng_1
82 : REAL(SP) :: r
83 :
84 : integer, parameter :: i8 = selected_int_kind(13)
85 : logical :: bigendian
86 :
87 123930000 : bigendian = (transfer(1_i8, 1) == 0)
88 :
89 : ! Return the next random numbers
90 :
91 : ! Marsaglia CONG algorithm
92 : if (bigendian) then
93 : ! Get low bytes from big-endian machine without overflow.
94 : s%seed = transfer(ishft(69069_i8*s%seed+1234567, bit_size(1)), 1)
95 : else
96 : ! Get low bytes from little-endian machine without overflow.
97 123930000 : s%seed = transfer(69069_i8*s%seed+1234567, 1)
98 : end if
99 :
100 : ! mod 32 bit overflow
101 123930000 : s%seed=mod(s%seed,2_ki9**30_ki9)
102 123930000 : r = s%seed*0.931322574615479E-09
103 :
104 : ! convert to range 0-1 (32 bit only)
105 : ! DJS2016: What is being done here is an intentional integer overflow and a test to
106 : ! see if this occured. Some compilers check for integer overflows during
107 : ! compilation (ie. gfortan), while others do not (ie. pgi and ifort). When
108 : ! using gfortran, you cannot use the overflow and test for overflow method,
109 : ! so we use sizeof(someInt) to determine wheter it is on 32 bit.
110 : !if ( i2_16*i2_16 .le. huge32 ) then
111 : if (digits(testInt) .le. 31) then
112 : !if (sizeof(testInt) .eq. 4) then
113 123930000 : r=r+1
114 123930000 : r=r-int(r)
115 : endif
116 123930000 : get_rng_1 = REAL(r, KIND = WP)
117 :
118 123930000 : END FUNCTION get_rng_1
119 :
120 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
121 : ! Create single random number N times
122 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
123 0 : FUNCTION get_rng_n(s, n) RESULT (r)
124 : integer,intent(inout) :: n
125 : TYPE(rng_state),INTENT(INOUT) :: s
126 : ! Return the next N random numbers
127 : REAL(WP), DIMENSION (n) :: r
128 :
129 : INTEGER :: i
130 :
131 0 : DO i = 1, N
132 0 : r(i) = get_rng_1(s)
133 : END DO
134 0 : END FUNCTION get_rng_n
135 :
136 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
137 : ! Create a vector of random numbers from a vector of input seeds
138 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
139 7894800 : FUNCTION get_rng_v(s) RESULT (r)
140 : ! Return the next N random numbers
141 : TYPE(rng_state), DIMENSION(:), INTENT(INOUT) :: s
142 : REAL(WP), DIMENSION (SIZE(S)) :: r
143 :
144 : INTEGER :: i
145 :
146 131824800 : DO i = 1, size(s)
147 131824800 : r(i) = get_rng_1(s(i))
148 : END DO
149 7894800 : END FUNCTION get_rng_v
150 :
151 0 : END MODULE mod_rng
|
