1 | MODULE storng |
---|
2 | !!====================================================================== |
---|
3 | !! *** MODULE storng *** |
---|
4 | !! Random number generator, used in NEMO stochastic parameterization |
---|
5 | !! |
---|
6 | !!===================================================================== |
---|
7 | !! History : 3.3 ! 2011-10 (J.-M. Brankart) Original code |
---|
8 | !!---------------------------------------------------------------------- |
---|
9 | |
---|
10 | !!---------------------------------------------------------------------- |
---|
11 | !! The module is based on (and includes) the |
---|
12 | !! 64-bit KISS (Keep It Simple Stupid) random number generator |
---|
13 | !! distributed by George Marsaglia : |
---|
14 | !! http://groups.google.com/group/comp.lang.fortran/ |
---|
15 | !! browse_thread/thread/a85bf5f2a97f5a55 |
---|
16 | !!---------------------------------------------------------------------- |
---|
17 | |
---|
18 | !!---------------------------------------------------------------------- |
---|
19 | !! kiss : 64-bit KISS random number generator (period ~ 2^250) |
---|
20 | !! kiss_seed : Define seeds for KISS random number generator |
---|
21 | !! kiss_state : Get current state of KISS random number generator |
---|
22 | !! kiss_save : Save current state of KISS (for future restart) |
---|
23 | !! kiss_load : Load the saved state of KISS |
---|
24 | !! kiss_reset : Reset the default seeds |
---|
25 | !! kiss_check : Check the KISS pseudo-random sequence |
---|
26 | !! kiss_uniform : Real random numbers with uniform distribution in [0,1] |
---|
27 | !! kiss_gaussian : Real random numbers with Gaussian distribution N(0,1) |
---|
28 | !! kiss_gamma : Real random numbers with Gamma distribution Gamma(k,1) |
---|
29 | !! kiss_sample : Select a random sample from a set of integers |
---|
30 | !! |
---|
31 | !! ---CURRENTLY NOT USED IN NEMO : |
---|
32 | !! kiss_save, kiss_load, kiss_check, kiss_gamma, kiss_sample |
---|
33 | !!---------------------------------------------------------------------- |
---|
34 | IMPLICIT NONE |
---|
35 | PRIVATE |
---|
36 | |
---|
37 | ! Public functions/subroutines |
---|
38 | PUBLIC :: kiss, kiss_seed, kiss_state, kiss_save, kiss_load, kiss_reset, kiss_check |
---|
39 | PUBLIC :: kiss_uniform, kiss_gaussian, kiss_gamma, kiss_sample |
---|
40 | |
---|
41 | ! Default/initial seeds |
---|
42 | INTEGER(KIND=8) :: x=1234567890987654321_8 |
---|
43 | INTEGER(KIND=8) :: y=362436362436362436_8 |
---|
44 | INTEGER(KIND=8) :: z=1066149217761810_8 |
---|
45 | INTEGER(KIND=8) :: w=123456123456123456_8 |
---|
46 | |
---|
47 | ! Parameters to generate real random variates |
---|
48 | REAL(KIND=8), PARAMETER :: huge64=9223372036854775808.0_8 ! +1 |
---|
49 | REAL(KIND=8), PARAMETER :: zero=0.0_8, half=0.5_8, one=1.0_8, two=2.0_8 |
---|
50 | |
---|
51 | ! Variables to store 2 Gaussian random numbers with current index (ig) |
---|
52 | INTEGER(KIND=8), SAVE :: ig=1 |
---|
53 | REAL(KIND=8), SAVE :: gran1, gran2 |
---|
54 | |
---|
55 | !!---------------------------------------------------------------------- |
---|
56 | !! NEMO/OPA 3.3 , NEMO Consortium (2010) |
---|
57 | !! $Id$ |
---|
58 | !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) |
---|
59 | !!---------------------------------------------------------------------- |
---|
60 | CONTAINS |
---|
61 | |
---|
62 | FUNCTION kiss() |
---|
63 | !! -------------------------------------------------------------------- |
---|
64 | !! *** FUNCTION kiss *** |
---|
65 | !! |
---|
66 | !! ** Purpose : 64-bit KISS random number generator |
---|
67 | !! |
---|
68 | !! ** Method : combine several random number generators: |
---|
69 | !! (1) Xorshift (XSH), period 2^64-1, |
---|
70 | !! (2) Multiply-with-carry (MWC), period (2^121+2^63-1) |
---|
71 | !! (3) Congruential generator (CNG), period 2^64. |
---|
72 | !! |
---|
73 | !! overall period: |
---|
74 | !! (2^250+2^192+2^64-2^186-2^129)/6 |
---|
75 | !! ~= 2^(247.42) or 10^(74.48) |
---|
76 | !! |
---|
77 | !! set your own seeds with 'kiss_seed' |
---|
78 | ! -------------------------------------------------------------------- |
---|
79 | IMPLICIT NONE |
---|
80 | INTEGER(KIND=8) :: kiss, t |
---|
81 | |
---|
82 | t = ISHFT(x,58) + w |
---|
83 | IF (s(x).eq.s(t)) THEN |
---|
84 | w = ISHFT(x,-6) + s(x) |
---|
85 | ELSE |
---|
86 | w = ISHFT(x,-6) + 1 - s(x+t) |
---|
87 | ENDIF |
---|
88 | x = t + x |
---|
89 | y = m( m( m(y,13_8), -17_8 ), 43_8 ) |
---|
90 | z = 6906969069_8 * z + 1234567_8 |
---|
91 | |
---|
92 | kiss = x + y + z |
---|
93 | |
---|
94 | CONTAINS |
---|
95 | |
---|
96 | FUNCTION s(k) |
---|
97 | INTEGER(KIND=8) :: s, k |
---|
98 | s = ISHFT(k,-63) |
---|
99 | END FUNCTION s |
---|
100 | |
---|
101 | FUNCTION m(k, n) |
---|
102 | INTEGER(KIND=8) :: m, k, n |
---|
103 | m = IEOR(k, ISHFT(k, n) ) |
---|
104 | END FUNCTION m |
---|
105 | |
---|
106 | END FUNCTION kiss |
---|
107 | |
---|
108 | |
---|
109 | SUBROUTINE kiss_seed(ix, iy, iz, iw) |
---|
110 | !! -------------------------------------------------------------------- |
---|
111 | !! *** ROUTINE kiss_seed *** |
---|
112 | !! |
---|
113 | !! ** Purpose : Define seeds for KISS random number generator |
---|
114 | !! |
---|
115 | !! -------------------------------------------------------------------- |
---|
116 | IMPLICIT NONE |
---|
117 | INTEGER(KIND=8) :: ix, iy, iz, iw |
---|
118 | |
---|
119 | x = ix |
---|
120 | y = iy |
---|
121 | z = iz |
---|
122 | w = iw |
---|
123 | |
---|
124 | END SUBROUTINE kiss_seed |
---|
125 | |
---|
126 | |
---|
127 | SUBROUTINE kiss_state(ix, iy, iz, iw) |
---|
128 | !! -------------------------------------------------------------------- |
---|
129 | !! *** ROUTINE kiss_state *** |
---|
130 | !! |
---|
131 | !! ** Purpose : Get current state of KISS random number generator |
---|
132 | !! |
---|
133 | !! -------------------------------------------------------------------- |
---|
134 | IMPLICIT NONE |
---|
135 | INTEGER(KIND=8) :: ix, iy, iz, iw |
---|
136 | |
---|
137 | ix = x |
---|
138 | iy = y |
---|
139 | iz = z |
---|
140 | iw = w |
---|
141 | |
---|
142 | END SUBROUTINE kiss_state |
---|
143 | |
---|
144 | |
---|
145 | SUBROUTINE kiss_reset() |
---|
146 | !! -------------------------------------------------------------------- |
---|
147 | !! *** ROUTINE kiss_reset *** |
---|
148 | !! |
---|
149 | !! ** Purpose : Reset the default seeds for KISS random number generator |
---|
150 | !! |
---|
151 | !! -------------------------------------------------------------------- |
---|
152 | IMPLICIT NONE |
---|
153 | |
---|
154 | x=1234567890987654321_8 |
---|
155 | y=362436362436362436_8 |
---|
156 | z=1066149217761810_8 |
---|
157 | w=123456123456123456_8 |
---|
158 | |
---|
159 | END SUBROUTINE kiss_reset |
---|
160 | |
---|
161 | |
---|
162 | SUBROUTINE kiss_check(check_type) |
---|
163 | !! -------------------------------------------------------------------- |
---|
164 | !! *** ROUTINE kiss_check *** |
---|
165 | !! |
---|
166 | !! ** Purpose : Check the KISS pseudo-random sequence |
---|
167 | !! |
---|
168 | !! ** Method : Check that it reproduces the correct sequence |
---|
169 | !! from the default seed |
---|
170 | !! |
---|
171 | !! -------------------------------------------------------------------- |
---|
172 | IMPLICIT NONE |
---|
173 | INTEGER(KIND=8) :: iter, niter, correct, iran |
---|
174 | CHARACTER(LEN=*) :: check_type |
---|
175 | LOGICAL :: print_success |
---|
176 | |
---|
177 | ! Save current state of KISS |
---|
178 | CALL kiss_save() |
---|
179 | ! Reset the default seed |
---|
180 | CALL kiss_reset() |
---|
181 | |
---|
182 | ! Select check type |
---|
183 | SELECT CASE(check_type) |
---|
184 | CASE('short') |
---|
185 | niter = 5_8 |
---|
186 | correct = 542381058189297533_8 |
---|
187 | print_success = .FALSE. |
---|
188 | CASE('long') |
---|
189 | niter = 100000000_8 |
---|
190 | correct = 1666297717051644203_8 ! Check provided by G. Marsaglia |
---|
191 | print_success = .TRUE. |
---|
192 | CASE('default') |
---|
193 | CASE DEFAULT |
---|
194 | STOP 'Bad check type in kiss_check' |
---|
195 | END SELECT |
---|
196 | |
---|
197 | ! Run kiss for the required number of iterations (niter) |
---|
198 | DO iter=1,niter |
---|
199 | iran = kiss() |
---|
200 | ENDDO |
---|
201 | |
---|
202 | ! Check that last iterate is correct |
---|
203 | IF (iran.NE.correct) THEN |
---|
204 | STOP 'Check failed: KISS internal error !!' |
---|
205 | ELSE |
---|
206 | IF (print_success) PRINT *, 'Check successful: 100 million calls to KISS OK' |
---|
207 | ENDIF |
---|
208 | |
---|
209 | ! Reload the previous state of KISS |
---|
210 | CALL kiss_load() |
---|
211 | |
---|
212 | END SUBROUTINE kiss_check |
---|
213 | |
---|
214 | |
---|
215 | SUBROUTINE kiss_save |
---|
216 | !! -------------------------------------------------------------------- |
---|
217 | !! *** ROUTINE kiss_save *** |
---|
218 | !! |
---|
219 | !! ** Purpose : Save current state of KISS random number generator |
---|
220 | !! |
---|
221 | !! -------------------------------------------------------------------- |
---|
222 | IMPLICIT NONE |
---|
223 | |
---|
224 | OPEN(UNIT=30,FILE='.kiss_restart') |
---|
225 | WRITE(30,*) x |
---|
226 | WRITE(30,*) y |
---|
227 | WRITE(30,*) z |
---|
228 | WRITE(30,*) w |
---|
229 | CLOSE(30) |
---|
230 | |
---|
231 | END SUBROUTINE kiss_save |
---|
232 | |
---|
233 | |
---|
234 | SUBROUTINE kiss_load |
---|
235 | !! -------------------------------------------------------------------- |
---|
236 | !! *** ROUTINE kiss_load *** |
---|
237 | !! |
---|
238 | !! ** Purpose : Load the saved state of KISS random number generator |
---|
239 | !! |
---|
240 | !! -------------------------------------------------------------------- |
---|
241 | IMPLICIT NONE |
---|
242 | LOGICAL :: filexists |
---|
243 | |
---|
244 | INQUIRE(FILE='.kiss_restart',EXIST=filexists) |
---|
245 | IF (filexists) THEN |
---|
246 | OPEN(UNIT=30,FILE='.kiss_restart') |
---|
247 | READ(30,*) x |
---|
248 | READ(30,*) y |
---|
249 | READ(30,*) z |
---|
250 | READ(30,*) w |
---|
251 | CLOSE(30) |
---|
252 | ENDIF |
---|
253 | |
---|
254 | END SUBROUTINE kiss_load |
---|
255 | |
---|
256 | |
---|
257 | SUBROUTINE kiss_uniform(uran) |
---|
258 | !! -------------------------------------------------------------------- |
---|
259 | !! *** ROUTINE kiss_uniform *** |
---|
260 | !! |
---|
261 | !! ** Purpose : Real random numbers with uniform distribution in [0,1] |
---|
262 | !! |
---|
263 | !! -------------------------------------------------------------------- |
---|
264 | IMPLICIT NONE |
---|
265 | REAL(KIND=8) :: uran |
---|
266 | |
---|
267 | uran = half * ( one + REAL(kiss(),8) / huge64 ) |
---|
268 | |
---|
269 | END SUBROUTINE kiss_uniform |
---|
270 | |
---|
271 | |
---|
272 | SUBROUTINE kiss_gaussian(gran) |
---|
273 | !! -------------------------------------------------------------------- |
---|
274 | !! *** ROUTINE kiss_gaussian *** |
---|
275 | !! |
---|
276 | !! ** Purpose : Real random numbers with Gaussian distribution N(0,1) |
---|
277 | !! |
---|
278 | !! ** Method : Generate 2 new Gaussian draws (gran1 and gran2) |
---|
279 | !! from 2 uniform draws on [-1,1] (u1 and u2), |
---|
280 | !! using the Marsaglia polar method |
---|
281 | !! (see Devroye, Non-Uniform Random Variate Generation, p. 235-236) |
---|
282 | !! |
---|
283 | !! -------------------------------------------------------------------- |
---|
284 | IMPLICIT NONE |
---|
285 | REAL(KIND=8) :: gran, u1, u2, rsq, fac |
---|
286 | |
---|
287 | IF (ig.EQ.1) THEN |
---|
288 | rsq = two |
---|
289 | DO WHILE ( (rsq.GE.one).OR. (rsq.EQ.zero) ) |
---|
290 | u1 = REAL(kiss(),8) / huge64 |
---|
291 | u2 = REAL(kiss(),8) / huge64 |
---|
292 | rsq = u1*u1 + u2*u2 |
---|
293 | ENDDO |
---|
294 | fac = SQRT(-two*LOG(rsq)/rsq) |
---|
295 | gran1 = u1 * fac |
---|
296 | gran2 = u2 * fac |
---|
297 | ENDIF |
---|
298 | |
---|
299 | ! Output one of the 2 draws |
---|
300 | IF (ig.EQ.1) THEN |
---|
301 | gran = gran1 ; ig = 2 |
---|
302 | ELSE |
---|
303 | gran = gran2 ; ig = 1 |
---|
304 | ENDIF |
---|
305 | |
---|
306 | END SUBROUTINE kiss_gaussian |
---|
307 | |
---|
308 | |
---|
309 | SUBROUTINE kiss_gamma(gamr,k) |
---|
310 | !! -------------------------------------------------------------------- |
---|
311 | !! *** ROUTINE kiss_gamma *** |
---|
312 | !! |
---|
313 | !! ** Purpose : Real random numbers with Gamma distribution Gamma(k,1) |
---|
314 | !! |
---|
315 | !! -------------------------------------------------------------------- |
---|
316 | IMPLICIT NONE |
---|
317 | REAL(KIND=8), PARAMETER :: p1 = 4.5_8 |
---|
318 | REAL(KIND=8), PARAMETER :: p2 = 2.50407739677627_8 ! 1+LOG(9/2) |
---|
319 | REAL(KIND=8), PARAMETER :: p3 = 1.38629436111989_8 ! LOG(4) |
---|
320 | REAL(KIND=8) :: gamr, k, u1, u2, b, c, d, xx, yy, zz, rr, ee |
---|
321 | LOGICAL :: accepted |
---|
322 | |
---|
323 | IF (k.GT.one) THEN |
---|
324 | ! Cheng's rejection algorithm |
---|
325 | ! (see Devroye, Non-Uniform Random Variate Generation, p. 413) |
---|
326 | b = k - p3 ; d = SQRT(two*k-one) ; c = k + d |
---|
327 | |
---|
328 | accepted=.FALSE. |
---|
329 | DO WHILE (.NOT.accepted) |
---|
330 | CALL kiss_uniform(u1) |
---|
331 | yy = LOG(u1/(one-u1)) / d ! Mistake in Devroye: "* k" instead of "/ d" |
---|
332 | xx = k * EXP(yy) |
---|
333 | rr = b + c * yy - xx |
---|
334 | CALL kiss_uniform(u2) |
---|
335 | zz = u1 * u1 * u2 |
---|
336 | |
---|
337 | accepted = rr .GE. (zz*p1-p2) |
---|
338 | IF (.NOT.accepted) accepted = rr .GE. LOG(zz) |
---|
339 | ENDDO |
---|
340 | |
---|
341 | gamr = xx |
---|
342 | |
---|
343 | ELSEIF (k.LT.one) THEN |
---|
344 | ! Rejection from the Weibull density |
---|
345 | ! (see Devroye, Non-Uniform Random Variate Generation, p. 415) |
---|
346 | c = one/k ; d = (one-k) * EXP( (k/(one-k)) * LOG(k) ) |
---|
347 | |
---|
348 | accepted=.FALSE. |
---|
349 | DO WHILE (.NOT.accepted) |
---|
350 | CALL kiss_uniform(u1) |
---|
351 | zz = -LOG(u1) |
---|
352 | xx = EXP( c * LOG(zz) ) |
---|
353 | CALL kiss_uniform(u2) |
---|
354 | ee = -LOG(u2) |
---|
355 | |
---|
356 | accepted = (zz+ee) .GE. (d+xx) ! Mistake in Devroye: "LE" instead of "GE" |
---|
357 | ENDDO |
---|
358 | |
---|
359 | gamr = xx |
---|
360 | |
---|
361 | ELSE |
---|
362 | ! Exponential distribution |
---|
363 | CALL kiss_uniform(u1) |
---|
364 | gamr = -LOG(u1) |
---|
365 | |
---|
366 | ENDIF |
---|
367 | |
---|
368 | END SUBROUTINE kiss_gamma |
---|
369 | |
---|
370 | |
---|
371 | SUBROUTINE kiss_sample(a,n,k) |
---|
372 | !! -------------------------------------------------------------------- |
---|
373 | !! *** ROUTINE kiss_sample *** |
---|
374 | !! |
---|
375 | !! ** Purpose : Select a random sample of size k from a set of n integers |
---|
376 | !! |
---|
377 | !! ** Method : The sample is output in the first k elements of a |
---|
378 | !! Set k equal to n to obtain a random permutation |
---|
379 | !! of the whole set of integers |
---|
380 | !! |
---|
381 | !! -------------------------------------------------------------------- |
---|
382 | IMPLICIT NONE |
---|
383 | INTEGER(KIND=8), DIMENSION(:) :: a |
---|
384 | INTEGER(KIND=8) :: n, k, i, j, atmp |
---|
385 | REAL(KIND=8) :: uran |
---|
386 | |
---|
387 | ! Select the sample using the swapping method |
---|
388 | ! (see Devroye, Non-Uniform Random Variate Generation, p. 612) |
---|
389 | DO i=1,k |
---|
390 | ! Randomly select the swapping element between i and n (inclusive) |
---|
391 | CALL kiss_uniform(uran) |
---|
392 | j = i - 1 + CEILING( REAL(n-i+1,8) * uran ) |
---|
393 | ! Swap elements i and j |
---|
394 | atmp = a(i) ; a(i) = a(j) ; a(j) = atmp |
---|
395 | ENDDO |
---|
396 | |
---|
397 | END SUBROUTINE kiss_sample |
---|
398 | |
---|
399 | END MODULE storng |
---|