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
|
