Branch data Line data Source code
1 : :
2 : : #if defined HAVE_CONFIG_H
3 : : #include "config.h"
4 : : #endif
5 : :
6 : : #include "abi_common.h"
7 : :
8 : : #if defined HAVE_FC_SHIFTLR
9 : : #define SHIFTL_ shiftl
10 : : #define SHIFTR_ shiftr
11 : : #else
12 : : #define SHIFTL_ lshift
13 : : #define SHIFTR_ rshift
14 : : #endif
15 : :
16 : : ! xoroshiro128plus method random number generator
17 : : ! adapted by hexu for usage in Abinit (downgrade to Fortran 90 and added
18 : : ! some functions )
19 : : ! (Fortran 2003 re-enabled)
20 : :
21 : : !** License for the xoroshiro128plus random number generator **
22 : : !
23 : : !Written in 2016 by David Blackman and Sebastiano Vigna (vigna@acm.org)
24 : : !Translated to Fortran 2008 by Jannis Teunissen
25 : : !
26 : : !To the extent possible under law, the author has dedicated all copyright
27 : : !and related and neighboring rights to this software to the public domain
28 : : !worldwide. This software is distributed without any warranty.
29 : : !
30 : : !See <http://creativecommons.org/publicdomain/zero/1.0/>.
31 : : !
32 : : ! << This is the original documentation of Blackman and Vigna
33 : : !This is the successor to xorshift128+. It is the fastest full-period
34 : : !generator passing BigCrush without systematic failures, but due to the
35 : : !relatively short period it is acceptable only for applications with a
36 : : !mild amount of parallelism; otherwise, use a xorshift1024* generator.
37 : : !
38 : : !Beside passing BigCrush, this generator passes the PractRand test suite
39 : : !up to (and included) 16TB, with the exception of binary rank tests,
40 : : !which fail due to the lowest bit being an LFSR; all other bits pass all
41 : : !tests. We suggest to use a sign test to extract a random Boolean value.
42 : : !
43 : : !Note that the generator uses a simulated rotate operation, which most C
44 : : !compilers will turn into a single instruction. In Java, you can use
45 : : !Long.rotateLeft(). In languages that do not make low-level rotation
46 : : !instructions accessible xorshift128+ could be faster.
47 : : !
48 : : !The state must be seeded so that it is not everywhere zero. If you have
49 : : !a 64-bit seed, we suggest to seed a splitmix64 generator and use its
50 : : !output to fill s.
51 : : ! =================
52 : : !
53 : :
54 : : !> Module for pseudo random number generation. The internal pseudo random
55 : : !> generator is the xoroshiro128plus method.
56 : :
57 : : module m_random_xoroshiro128plus
58 : :
59 : : use defs_basis
60 : : use m_profiling_abi
61 : :
62 : : implicit none
63 : : private
64 : :
65 : : ! A 64 bit floating point type
66 : : !integer, parameter :: dp = kind(0.0d0)
67 : :
68 : : ! A 32 bit integer type
69 : : integer, parameter :: i4 = selected_int_kind(9)
70 : :
71 : : ! A 64 bit integer type
72 : : integer, parameter :: i8 = selected_int_kind(18)
73 : :
74 : : !> Random number generator type, which contains the state
75 : : type rng_t
76 : : !> The rng state (always use your own seed)
77 : : integer(i8), private :: s(2) = [123456789_i8, 987654321_i8]
78 : : integer(i8), private :: separator(32) ! Separate cache lines (parallel use)
79 : :
80 : : real(dp) :: residual =0.0d0 ! for saving residual in normal function.
81 : : logical :: has_residual = .False.
82 : : contains
83 : : procedure, non_overridable :: set_seed ! Seed the generator
84 : : procedure, non_overridable :: jump ! Jump function (see below)
85 : : procedure, non_overridable :: rand_int4 ! 4-byte random integer
86 : : procedure, non_overridable :: rand_int8 ! 8-byte random integer
87 : : procedure, non_overridable :: rand_unif_01 ! Uniform (0,1] real
88 : : procedure, non_overridable :: rand_unif_01_array ! Uniform (0,1] real
89 : : procedure, non_overridable :: rand_two_normals ! Two normal(0,1) samples
90 : : procedure, non_overridable :: rand_normal ! Two normal(0,1) samples
91 : : procedure, non_overridable :: rand_normal_array ! Two normal(0,1) samples
92 : : procedure, non_overridable :: rand_poisson ! Sample from Poisson-dist.
93 : : procedure, non_overridable :: rand_circle ! Sample on a rand_circle
94 : : procedure, non_overridable :: rand_sphere ! Sample on a rand_sphere
95 : : procedure, non_overridable :: rand_choice ! select from 1 to N randomly
96 : : procedure, non_overridable :: next ! Internal method
97 : : end type rng_t
98 : :
99 : : !> Parallel random number generator type
100 : : type prng_t
101 : : type(rng_t), allocatable :: rngs(:)
102 : : contains
103 : : procedure, non_overridable :: initialize => init_parallel
104 : : procedure, non_overridable :: finalize => finalize_parallel
105 : : end type prng_t
106 : :
107 : :
108 : : public :: rng_t
109 : : public :: prng_t
110 : : public :: set_seed
111 : : public :: jump
112 : : public :: rand_int4
113 : : public :: rand_int8
114 : : public :: rand_unif_01
115 : : public :: rand_two_normals
116 : : public :: rand_normal
117 : : public :: rand_normal_array
118 : : public :: rand_poisson
119 : : public :: rand_circle
120 : : public :: rand_sphere
121 : :
122 : :
123 : : contains
124 : :
125 : : !> Initialize a collection of rng's for parallel use
126 : : !> this can be used with openmp, each thread need an
127 : : !> different rng.
128 : 0 : subroutine init_parallel(self, n_proc, seed)
129 : :
130 : : class(prng_t), intent(inout) :: self
131 : : integer(i8), optional, intent(in) :: seed(2)
132 : : integer, intent(in) :: n_proc
133 : : integer :: n
134 : :
135 [ # # # # : 0 : ABI_MALLOC(self%rngs, (n_proc))
# # # # #
# # # # #
# # # # ]
136 [ # # ]: 0 : if (present(seed)) then
137 : 0 : call self%rngs(1)%set_seed(seed)
138 : : endif
139 [ # # ]: 0 : do n = 2, n_proc
140 : 0 : self%rngs(n) = self%rngs(n-1)
141 : 0 : call self%rngs(n)%jump()
142 : : end do
143 : 0 : end subroutine init_parallel
144 : :
145 : :
146 : : !> Finalize prng_t
147 : 0 : subroutine finalize_parallel(self)
148 : : class(prng_t), intent(inout) :: self
149 [ # # ]: 0 : ABI_FREE(self%rngs)
150 : 0 : end subroutine finalize_parallel
151 : :
152 : : !> Set a seed for the rng
153 : 6 : subroutine set_seed(self, the_seed)
154 : :
155 : : class(rng_t), intent(inout) :: self
156 : : integer(i8), intent(in) :: the_seed(2)
157 : :
158 [ + + ]: 18 : self%s = the_seed
159 : :
160 : : ! Simulate calls to next() to improve randomness of first number
161 : : !call self%jump()
162 : 6 : call jump(self)
163 : 6 : end subroutine set_seed
164 : :
165 : : ! This is the jump function for the generator. It is equivalent
166 : : ! to 2^64 calls to next(); it can be used to generate 2^64
167 : : ! non-overlapping subsequences for parallel computations.
168 : 6 : subroutine jump(self)
169 : :
170 : : class(rng_t), intent(inout) :: self
171 : : integer :: i, b
172 : : integer(i8) :: t(2), dummy
173 : :
174 : : ! The signed equivalent of the unsigned constants
175 : : integer(i8), parameter :: jmp_c(2) = &
176 : : (/-4707382666127344949_i8, -2852180941702784734_i8/)
177 : :
178 : 6 : t = 0
179 [ + + ]: 18 : do i = 1, 2
180 [ + + ]: 786 : do b = 0, 63
181 [ + + ]: 768 : if (iand(jmp_c(i), SHIFTL_(1_i8, b)) /= 0) then
182 [ + + ]: 1152 : t = ieor(t, self%s)
183 : : end if
184 : : !dummy = self%next()
185 : 780 : dummy = next(self)
186 : : end do
187 : : end do
188 : :
189 [ + + ]: 18 : self%s = t
190 : 6 : end subroutine jump
191 : :
192 : : !> Return 4-byte integer
193 : 0 : integer(i4) function rand_int4(self)
194 : :
195 : : class(rng_t), intent(inout) :: self
196 : : !rand_int4 = int(self%next(), i4)
197 : 0 : rand_int4 = int(next(self), i4)
198 : 0 : end function rand_int4
199 : :
200 : : !> Return 8-byte integer
201 : 0 : integer(i8) function rand_int8(self)
202 : :
203 : : class(rng_t), intent(inout) :: self
204 : : !rand_int8 = self%next()
205 : 0 : rand_int8 = next(self)
206 : 0 : end function rand_int8
207 : :
208 : : !> Get a uniform [0,1) random real (double precision)
209 : 1697616 : real(dp) function rand_unif_01(self)
210 : :
211 : : class(rng_t), intent(inout) :: self
212 : : integer(i8) :: x
213 : : real(dp) :: tmp
214 : :
215 : : !x = self%next()
216 : 848808 : x = next(self)
217 : 1697616 : x = ior(SHIFTL_(1023_i8, 52), SHIFTR_(x, 12))
218 : 1697616 : rand_unif_01 = transfer(x, tmp) - 1.0_dp
219 : 0 : end function rand_unif_01
220 : :
221 : :
222 : 4 : subroutine rand_unif_01_array(self, output, size_array)
223 : :
224 : : class(rng_t), intent(inout) :: self
225 : : integer, intent(in) :: size_array
226 : : real(dp), intent(inout) :: output(size_array)
227 : : integer(i8) :: x, i
228 : : real(dp) :: tmp
229 : :
230 : : !x = self%next()
231 [ + + ]: 3348 : do i=1, size_array
232 : 3344 : x = next(self)
233 : 3344 : x = ior(SHIFTL_(1023_i8, 52), SHIFTR_(x, 12))
234 : 3348 : output(i) = transfer(x, tmp) - 1.0_dp
235 : : enddo
236 : 4 : end subroutine rand_unif_01_array
237 : :
238 : : !> Return two normal random variates with mean 0 and variance 1.
239 : : !> http://en.wikipedia.org/wiki/Marsaglia_polar_method
240 [ # # ]: 0 : function rand_two_normals(self) result(rands)
241 : :
242 : : class(rng_t), intent(inout) :: self
243 : : real(dp) :: rands(2), sum_sq
244 : : do
245 : : !rands(1) = 2 * self%rand_unif_01() - 1
246 : 0 : rands(1) = 2 * rand_unif_01(self) - 1
247 : : !rands(2) = 2 * self%rand_unif_01() - 1
248 : 0 : rands(2) = 2 * rand_unif_01(self) - 1
249 [ # # ]: 0 : sum_sq = sum(rands**2)
250 [ # # # # ]: 0 : if (sum_sq < 1.0_dp .and. sum_sq > 0.0_dp) exit
251 : : end do
252 [ # # ]: 0 : rands = rands * sqrt(-2 * log(sum_sq) / sum_sq)
253 : 0 : end function rand_two_normals
254 : :
255 : : !> Return one normal random variates with mean 0 and variance 1.
256 : : !> http://en.wikipedia.org/wiki/Marsaglia_polar_method
257 : 1332464 : function rand_normal(self) result(r)
258 : :
259 : : class(rng_t), intent(inout) :: self
260 : : real(dp) :: r, rands(2), sum_sq
261 [ + + ]: 1332464 : if(self%has_residual) then
262 : 666232 : self%has_residual=.False.
263 : 666232 : r=self%residual
264 : : else
265 : : do
266 : : !rands(1) = 2 * self%rand_unif_01() - 1
267 : 848808 : rands(1) = 2 * rand_unif_01(self) - 1
268 : : !rands(2) = 2 * self%rand_unif_01() - 1
269 : 848808 : rands(2) = 2 * rand_unif_01(self) - 1
270 [ + + ]: 2546424 : sum_sq = sum(rands**2)
271 [ + + - + ]: 848808 : if (sum_sq < 1.0_dp .and. sum_sq > 0.0_dp) exit
272 : : end do
273 [ + + ]: 1998696 : rands = rands * sqrt(-2 * log(sum_sq) / sum_sq)
274 : 666232 : r=rands(1)
275 : 666232 : self%has_residual=.True.
276 : 666232 : self%residual=rands(2)
277 : : endif
278 : 1332464 : end function rand_normal
279 : :
280 : 2029 : subroutine rand_normal_array(self, a, asize)
281 : :
282 : : class(rng_t), intent(inout) :: self
283 : : integer ,intent(in) :: asize
284 : : real(dp), intent(inout) :: a(asize)
285 : : integer :: i
286 [ + + ]: 1334493 : do i=1, asize
287 : 1334493 : a(i)=rand_normal(self)
288 : : enddo
289 : 2029 : end subroutine rand_normal_array
290 : :
291 : :
292 : : !> Return Poisson random variate with rate lambda. Works well for lambda < 30
293 : : !> or so. For lambda >> 1 it can produce wrong results due to roundoff error.
294 : 0 : function rand_poisson(self, lambda) result(rr)
295 : :
296 : : class(rng_t), intent(inout) :: self
297 : : real(dp), intent(in) :: lambda
298 : : integer(i4) :: rr
299 : : real(dp) :: expl, p
300 : :
301 : 0 : expl = exp(-lambda)
302 : 0 : rr = 0
303 : 0 : p = self%rand_unif_01()
304 : :
305 [ # # ]: 0 : do while (p > expl)
306 : 0 : rr = rr + 1
307 : : !p = p * self%rand_unif_01()
308 : 0 : p = p * rand_unif_01(self)
309 : : end do
310 : 0 : end function rand_poisson
311 : :
312 : : !> Sample point on a rand_circle with given radius
313 [ # # ]: 0 : function rand_circle(self, radius) result(xy)
314 : :
315 : : class(rng_t), intent(inout) :: self
316 : : real(dp), intent(in) :: radius
317 : : real(dp) :: rands(2), xy(2)
318 : : real(dp) :: sum_sq
319 : :
320 : : ! Method for uniform sampling on rand_circle
321 : : do
322 : : !rands(1) = 2 * self%rand_unif_01() - 1
323 : 0 : rands(1) = 2 * rand_unif_01(self) - 1
324 : : !rands(2) = 2 * self%rand_unif_01() - 1
325 : 0 : rands(2) = 2 * rand_unif_01(self) - 1
326 [ # # ]: 0 : sum_sq = sum(rands**2)
327 [ # # ]: 0 : if (sum_sq <= 1) exit
328 : : end do
329 : :
330 : 0 : xy(1) = (rands(1)**2 - rands(2)**2) / sum_sq
331 : 0 : xy(2) = 2 * rands(1) * rands(2) / sum_sq
332 [ # # ]: 0 : xy = xy * radius
333 : 0 : end function rand_circle
334 : :
335 : : !> Sample point on a rand_sphere with given radius
336 [ # # ]: 0 : function rand_sphere(self, radius) result(xyz)
337 : :
338 : : class(rng_t), intent(inout) :: self
339 : : real(dp), intent(in) :: radius
340 : : real(dp) :: rands(2), xyz(3)
341 : : real(dp) :: sum_sq, tmp_sqrt
342 : :
343 : : ! Marsaglia method for uniform sampling on rand_sphere
344 : : do
345 : 0 : rands(1) = 2 * self%rand_unif_01() - 1
346 : 0 : rands(2) = 2 * self%rand_unif_01() - 1
347 [ # # ]: 0 : sum_sq = sum(rands**2)
348 [ # # ]: 0 : if (sum_sq <= 1) exit
349 : : end do
350 : :
351 : 0 : tmp_sqrt = sqrt(1 - sum_sq)
352 [ # # ]: 0 : xyz(1:2) = 2 * rands(1:2) * tmp_sqrt
353 : 0 : xyz(3) = 1 - 2 * sum_sq
354 [ # # ]: 0 : xyz = xyz * radius
355 : 0 : end function rand_sphere
356 : :
357 : 0 : function rand_choice(self, N) result(i)
358 : : class(rng_t), intent(inout) :: self
359 : : integer, intent(in) :: N
360 : : integer :: i
361 : 0 : i=int(self%rand_unif_01()*N)+1
362 : 0 : end function rand_choice
363 : :
364 : : !> Interal routine: get the next value (returned as 64 bit signed integer)
365 : 1701728 : function next(self) result(res)
366 : :
367 : : class(rng_t), intent(inout) :: self
368 : : integer(i8) :: res
369 : : integer(i8) :: t(2)
370 : :
371 [ + + ]: 5105184 : t = self%s
372 : 1701728 : res = t(1) + t(2)
373 : 1701728 : t(2) = ieor(t(1), t(2))
374 : 1701728 : self%s(1) = ieor(ieor(rotl(t(1), 55), t(2)), SHIFTL_(t(2), 14))
375 : 1701728 : self%s(2) = rotl(t(2), 36)
376 : 1701728 : end function next
377 : :
378 : : !> Helper function for next()
379 : 1701728 : pure function rotl(x, k) result(res)
380 : :
381 : : integer(i8), intent(in) :: x
382 : : integer, intent(in) :: k
383 : : integer(i8) :: res
384 : :
385 : 1701728 : res = ior(SHIFTL_(x, k), SHIFTR_(x, 64 - k))
386 : : end function rotl
387 : :
388 [ # # # # : 0 : end module m_random_xoroshiro128plus
# # # # #
# # # ]
|
