Home

Context Navigation

source: palm/trunk/SOURCE/random_generator_parallel_mod.f90 @ 2236

Last change on this file since 2236 was 2173, checked in by knoop, 8 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 18.4 KB

Line
1	!> @file random_generator_parallel_mod.f90
2	!------------------------------------------------------------------------------!
3	! This file is part of PALM.
4	!
5	! PALM is free software: you can redistribute it and/or modify it under the
6	! terms of the GNU General Public License as published by the Free Software
7	! Foundation, either version 3 of the License, or (at your option) any later
8	! version.
9	!
10	! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
11	! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
12	! A PARTICULAR PURPOSE. See the GNU General Public License for more details.
13	!
14	! You should have received a copy of the GNU General Public License along with
15	! PALM. If not, see <http://www.gnu.org/licenses/>.
16	!
17	! Copyright 1997-2017 Leibniz Universitaet Hannover
18	!------------------------------------------------------------------------------!
19	!
20	! Current revisions:
21	! -----------------
22	!
23	!
24	! Former revisions:
25	! -----------------
26	! $Id: random_generator_parallel_mod.f90 2173 2017-03-08 15:56:54Z suehring $
27	!
28	! 2172 2017-03-08 15:55:25Z knoop
29	! Bugfix: added global initialization routine and removed global array
30	!
31	! 2144 2017-02-07 15:32:45Z knoop
32	! Bugfix: compile time parameter replacement prevented
33	!
34	! 2000 2016-08-20 18:09:15Z knoop
35	! Forced header and separation lines into 80 columns
36	!
37	! 1850 2016-04-08 13:29:27Z maronga
38	! Module renamed
39	!
40	!
41	! 1682 2015-10-07 23:56:08Z knoop
42	! Code annotations made doxygen readable
43	!
44	! 1400 2014-05-09 14:03:54Z knoop
45	! Initial revision
46	!
47	!
48	! Description:
49	! ------------
50	!> This module contains and supports the random number generating routine ran_parallel.
51	!> ran_parallel returns a uniform random deviate between 0.0 and 1.0
52	!> (exclusive of the end point values).
53	!> Additionally it provides the generator with five integer for use as initial state space.
54	!> The first tree integers (iran, jran, kran) are maintained as non negative values,
55	!> while the last two (mran, nran) have 32-bit nonzero values.
56	!> Also provided by this module is support for initializing or reinitializing
57	!> the state space to a desired standard sequence number, hashing the initial
58	!> values to random values, and allocating and deallocating the internal workspace
59	!> Random number generator, produces numbers equally distributed in interval
60	!>
61	!> This routine is taken from the "numerical recipies vol. 2"
62	!------------------------------------------------------------------------------!
63	MODULE random_generator_parallel
64
65
66	USE kinds
67
68	IMPLICIT NONE
69
70	PRIVATE
71	PUBLIC random_number_parallel, random_seed_parallel, random_dummy, &
72	id_random_array, seq_random_array, init_parallel_random_generator
73
74	INTEGER(isp), SAVE :: lenran=0 !<
75	INTEGER(isp), SAVE :: seq=0 !<
76	INTEGER(isp), SAVE :: iran0 !<
77	INTEGER(isp), SAVE :: jran0 !<
78	INTEGER(isp), SAVE :: kran0 !<
79	INTEGER(isp), SAVE :: mran0 !<
80	INTEGER(isp), SAVE :: nran0 !<
81	INTEGER(isp), SAVE :: rans !<
82
83	INTEGER(isp), DIMENSION(:, :), POINTER, SAVE :: ranseeds !<
84
85	INTEGER(isp), DIMENSION(:), POINTER, SAVE :: iran !<
86	INTEGER(isp), DIMENSION(:), POINTER, SAVE :: jran !<
87	INTEGER(isp), DIMENSION(:), POINTER, SAVE :: kran !<
88	INTEGER(isp), DIMENSION(:), POINTER, SAVE :: mran !<
89	INTEGER(isp), DIMENSION(:), POINTER, SAVE :: nran !<
90	INTEGER(isp), DIMENSION(:), POINTER, SAVE :: ranv !<
91
92
93
94	INTEGER(isp), DIMENSION(:,:), ALLOCATABLE :: id_random_array !<
95	INTEGER(isp), DIMENSION(:,:,:), ALLOCATABLE :: seq_random_array !<
96
97	REAL(wp), SAVE :: amm !<
98
99	REAL(wp) :: random_dummy=0.0 !<
100
101	INTERFACE init_parallel_random_generator
102	MODULE PROCEDURE init_parallel_random_generator
103	END INTERFACE
104
105	INTERFACE random_number_parallel
106	MODULE PROCEDURE ran0_s
107	END INTERFACE
108
109	INTERFACE random_seed_parallel
110	MODULE PROCEDURE random_seed_parallel
111	END INTERFACE
112
113	INTERFACE ran_hash
114	MODULE PROCEDURE ran_hash_v
115	END INTERFACE
116
117	INTERFACE reallocate
118	MODULE PROCEDURE reallocate_iv,reallocate_im
119	END INTERFACE
120
121	INTERFACE arth
122	MODULE PROCEDURE arth_i
123	END INTERFACE
124
125	CONTAINS
126
127	!------------------------------------------------------------------------------!
128	! Description:
129	! ------------
130	!> Initialize the parallel random number generator for a specific subdomain
131	!------------------------------------------------------------------------------!
132	SUBROUTINE init_parallel_random_generator(nx, ny, nys, nyn, nxl, nxr)
133
134	USE kinds
135
136	USE control_parameters, &
137	ONLY: ensemble_member_nr
138
139	IMPLICIT NONE
140
141	INTEGER(isp), INTENT(IN) :: nx !<
142	INTEGER(isp), INTENT(IN) :: ny !<
143	INTEGER(isp), INTENT(IN) :: nys !<
144	INTEGER(isp), INTENT(IN) :: nyn !<
145	INTEGER(isp), INTENT(IN) :: nxl !<
146	INTEGER(isp), INTENT(IN) :: nxr !<
147
148	INTEGER(iwp) :: i !<
149	INTEGER(iwp) :: j !<
150
151	!-- Allocate ID-array and state-space-array
152	ALLOCATE ( seq_random_array(5,nys:nyn,nxl:nxr) )
153	ALLOCATE ( id_random_array(nys:nyn,nxl:nxr) )
154	seq_random_array = 0
155	id_random_array = 0
156
157	!-- Asigning an ID to every vertical gridpoint column
158	!-- dependig on the ensemble run number.
159	DO i=nxl, nxr
160	DO j=nys, nyn
161	id_random_array(j,i) = j(nx+1.0_wp) + i + 1.0_wp + 1E6 &
162	( ensemble_member_nr - 1000 )
163	ENDDO
164	ENDDO
165	!-- Initializing with random_seed_parallel for every vertical
166	!-- gridpoint column.
167	random_dummy=0
168	DO i = nxl, nxr
169	DO j = nys, nyn
170	CALL random_seed_parallel (random_sequence=id_random_array(j, i))
171	CALL random_number_parallel (random_dummy)
172	CALL random_seed_parallel (get=seq_random_array(:, j, i))
173	ENDDO
174	ENDDO
175
176	END SUBROUTINE init_parallel_random_generator
177
178	!------------------------------------------------------------------------------!
179	! Description:
180	! ------------
181	!> Lagged Fibonacci generator combined with a Marsaglia shiftsequence.
182	!> Returns as harvest a uniform random deviate between 0.0 and 1.0 (exclusive of the end point values).
183	!> This generator has the same calling and initialization conventions as Fortran 90's random_number routine.
184	!> Use random_seed_parallel to initialize or reinitialize to a particular sequence.
185	!> The period of this generator is about 2.0 x 10^28, and it fully vectorizes.
186	!> Validity of the integer model assumed by this generator is tested at initialization.
187	!------------------------------------------------------------------------------!
188	SUBROUTINE ran0_s(harvest)
189
190	USE kinds
191
192	IMPLICIT NONE
193
194	REAL(wp), INTENT(OUT) :: harvest !<
195
196	IF (lenran < 1) CALL ran_init(1) !- Initialization routine in ran_state.
197
198	!- Update Fibonacci generator, which has period p^2 + p + 1, p = 2^31 - 69.
199	rans = iran0 - kran0
200
201	IF (rans < 0) rans = rans + 2147483579_isp
202
203	iran0 = jran0
204	jran0 = kran0
205	kran0 = rans
206
207	nran0 = ieor( nran0, ishft (nran0, 13) ) !- Update Marsaglia shift sequence with period 2^32 - 1.
208	nran0 = ieor( nran0, ishft (nran0, -17) )
209	nran0 = ieor( nran0, ishft (nran0, 5) )
210
211	rans = ieor( nran0, rans ) !- Combine the generators.
212
213	harvest = amm * merge( rans, not(rans), rans < 0 ) !- Make the result positive definite (note that amm is negative).
214
215	END SUBROUTINE ran0_s
216
217	!------------------------------------------------------------------------------!
218	! Description:
219	! ------------
220	!> Initialize or reinitialize the random generator state space to vectors of size length.
221	!> The saved variable seq is hashed (via calls to the module routine ran_hash)
222	!> to create unique starting seeds, different for each vector component.
223	!------------------------------------------------------------------------------!
224	SUBROUTINE ran_init( length )
225
226	USE kinds
227
228	IMPLICIT NONE
229
230	INTEGER(isp), INTENT(IN) :: length !<
231
232	INTEGER(isp) :: hg !<
233	INTEGER(isp) :: hgm !<
234	INTEGER(isp) :: hgng !<
235
236	INTEGER(isp) :: new !<
237	INTEGER(isp) :: j !<
238	INTEGER(isp) :: hgt !<
239
240	IF ( length < lenran ) RETURN !- Simply return if enough space is already allocated.
241
242	hg=huge(1_isp)
243	hgm=-hg
244	hgng=hgm-1
245	hgt = hg
246
247	!- The following lines check that kind value isp is in fact a 32-bit integer
248	!- with the usual properties that we expect it to have (under negation and wrap-around addition).
249	!- If all of these tests are satisfied, then the routines that use this module are portable,
250	!- even though they go beyond Fortran 90's integer model.
251
252	IF ( hg /= 2147483647 ) CALL ran_error('ran_init: arith assump 1 fails')
253	IF ( hgng >= 0 ) CALL ran_error('ran_init: arith assump 2 fails')
254	IF ( hgt+1 /= hgng ) CALL ran_error('ran_init: arith assump 3 fails')
255	IF ( not(hg) >= 0 ) CALL ran_error('ran_init: arith assump 4 fails')
256	IF ( not(hgng) < 0 ) CALL ran_error('ran_init: arith assump 5 fails')
257	IF ( hg+hgng >= 0 ) CALL ran_error('ran_init: arith assump 6 fails')
258	IF ( not(-1_isp) < 0 ) CALL ran_error('ran_init: arith assump 7 fails')
259	IF ( not(0_isp) >= 0 ) CALL ran_error('ran_init: arith assump 8 fails')
260	IF ( not(1_isp) >= 0 ) CALL ran_error('ran_init: arith assump 9 fails')
261
262	IF ( lenran > 0) THEN !- Reallocate space, or ...
263
264	ranseeds => reallocate( ranseeds, length, 5)
265	ranv => reallocate( ranv, length-1)
266	new = lenran+1
267
268	ELSE !- allocate space.
269
270	ALLOCATE(ranseeds(length,5))
271	ALLOCATE(ranv(length-1))
272	new = 1 !- Index of first location not yet initialized.
273	amm = nearest(1.0_wp,-1.0_wp)/hgng
274	!- Use of nearest is to ensure that returned random deviates are strictly lessthan 1.0.
275	IF (ammhgng >= 1.0 .or. ammhgng <= 0.0) &
276	CALL ran_error('ran_init: arith assump 10 fails')
277
278	END IF
279
280	!- Set starting values, unique by seq and vector component.
281	ranseeds(new:,1) = seq
282	ranseeds(new:,2:5)=spread(arth(new,1,size(ranseeds(new:,1))),2,4)
283
284	DO j=1,4 !- Hash them.
285	CALL ran_hash(ranseeds(new:,j), ranseeds(new:,j+1))
286	END DO
287
288	WHERE (ranseeds (new: ,1:3) < 0) &
289	ranseeds(new: ,1:3)=not(ranseeds(new: ,1:3)) !- Enforce nonnegativity.
290
291	WHERE (ranseeds(new: ,4:5) == 0) ranseeds(new: ,4:5)=1 !- Enforce nonzero.
292
293	IF (new == 1) THEN !- Set scalar seeds.
294
295	iran0 = ranseeds(1,1)
296	jran0 = ranseeds(1,2)
297	kran0 = ranseeds(1,3)
298	mran0 = ranseeds(1,4)
299	nran0 = ranseeds(1,5)
300	rans = nran0
301
302	END IF
303
304	IF (length > 1) THEN !- Point to vector seeds.
305
306	iran => ranseeds(2:,1)
307	jran => ranseeds(2:,2)
308	kran => ranseeds(2:,3)
309	mran => ranseeds(2:,4)
310	nran => ranseeds(2:,5)
311	ranv = nran
312
313	END IF
314
315	lenran = length
316
317	END SUBROUTINE ran_init
318
319	!------------------------------------------------------------------------------!
320	! Description:
321	! ------------
322	!> User interface to release the workspace used by the random number routines.
323	!------------------------------------------------------------------------------!
324	SUBROUTINE ran_deallocate
325
326	IF ( lenran > 0 ) THEN
327
328	DEALLOCATE(ranseeds, ranv)
329	NULLIFY(ranseeds, ranv, iran, jran, kran, mran, nran)
330	lenran = 0
331
332	END IF
333
334	END SUBROUTINE ran_deallocate
335
336	!------------------------------------------------------------------------------!
337	! Description:
338	! ------------
339	!> User interface for seeding the random number routines.
340	!> Syntax is exactly like Fortran 90's random_seed routine,
341	!> with one additional argument keyword: random_sequence, set to any integer
342	!> value, causes an immediate new initialization, seeded by that integer.
343	!------------------------------------------------------------------------------!
344	SUBROUTINE random_seed_parallel( random_sequence, state_size, put, get )
345
346	IMPLICIT NONE
347
348	INTEGER(isp), OPTIONAL, INTENT(IN) :: random_sequence !<
349	INTEGER(isp), OPTIONAL, INTENT(OUT) :: state_size !<
350
351	INTEGER(isp), DIMENSION(:), OPTIONAL, INTENT(IN) :: put !<
352	INTEGER(isp), DIMENSION(:), OPTIONAL, INTENT(OUT) :: get !<
353
354	IF ( present(state_size) ) THEN
355
356	state_size = 5 * lenran
357
358	ELSE IF ( present(put) ) THEN
359
360	IF ( lenran == 0 ) RETURN
361
362	ranseeds = reshape( put,shape(ranseeds) )
363
364	WHERE (ranseeds(:,1:3) < 0) ranseeds(: ,1:3) = not(ranseeds(: ,1:3))
365	!- Enforce nonnegativity and nonzero conditions on any user-supplied seeds.
366
367	WHERE (ranseeds(:,4:5) == 0) ranseeds(:,4:5) = 1
368
369	iran0 = ranseeds(1,1)
370	jran0 = ranseeds(1,2)
371	kran0 = ranseeds(1,3)
372	mran0 = ranseeds(1,4)
373	nran0 = ranseeds(1,5)
374
375	ELSE IF ( present(get) ) THEN
376
377	IF (lenran == 0) RETURN
378
379	ranseeds(1,1:5) = (/ iran0,jran0,kran0,mran0,nran0 /)
380	get = reshape( ranseeds, shape(get) )
381
382	ELSE IF ( present(random_sequence) ) THEN
383
384	CALL ran_deallocate
385	seq = random_sequence
386
387	END IF
388
389	END SUBROUTINE random_seed_parallel
390
391	!------------------------------------------------------------------------------!
392	! Description:
393	! ------------
394	!> DES-like hashing of two 32-bit integers, using shifts,
395	!> xor's, and adds to make the internal nonlinear function.
396	!------------------------------------------------------------------------------!
397	SUBROUTINE ran_hash_v( il, ir )
398
399	IMPLICIT NONE
400
401	INTEGER(isp), DIMENSION(:), INTENT(INOUT) :: il !<
402	INTEGER(isp), DIMENSION(:), INTENT(INOUT) :: ir !<
403
404	INTEGER(isp), DIMENSION(size(il)) :: is !<
405
406	INTEGER(isp) :: j !<
407
408	DO j=1,4
409
410	is = ir
411	ir = ieor( ir, ishft(ir,5) ) + 1422217823
412	ir = ieor( ir, ishft(ir,-16) ) + 1842055030
413	ir = ieor( ir, ishft(ir,9) ) + 80567781
414	ir = ieor( il, ir )
415	il = is
416	END DO
417
418	END SUBROUTINE ran_hash_v
419
420	!------------------------------------------------------------------------------!
421	! Description:
422	! ------------
423	!> User interface to process error-messages
424	!> produced by the random_number_parallel module
425	!------------------------------------------------------------------------------!
426	SUBROUTINE ran_error(string)
427
428	CHARACTER(LEN=*), INTENT(IN) :: string !<
429
430	write (,) 'Error in module random_number_parallel: ',string
431
432	STOP 'Program terminated by ran_error'
433
434	END SUBROUTINE ran_error
435
436	!------------------------------------------------------------------------------!
437	! Description:
438	! ------------
439	!> Reallocates the generators state space "ranseeds" to vectors of size length.
440	!------------------------------------------------------------------------------!
441	FUNCTION reallocate_iv( p, n )
442
443	INTEGER(isp), DIMENSION(:), POINTER :: p !<
444	INTEGER(isp), DIMENSION(:), POINTER :: reallocate_iv !<
445
446	INTEGER(isp), INTENT(IN) :: n !<
447
448	INTEGER(isp) :: nold !<
449	INTEGER(isp) :: ierr !<
450
451	ALLOCATE(reallocate_iv(n),stat=ierr)
452
453	IF (ierr /= 0) CALL &
454	ran_error('reallocate_iv: problem in attempt to allocate memory')
455
456	IF (.not. associated(p)) RETURN
457
458	nold = size(p)
459
460	reallocate_iv(1:min(nold,n)) = p(1:min(nold,n))
461
462	DEALLOCATE(p)
463
464	END FUNCTION reallocate_iv
465
466	FUNCTION reallocate_im( p, n, m )
467
468	INTEGER(isp), DIMENSION(:,:), POINTER :: p !<
469	INTEGER(isp), DIMENSION(:,:), POINTER :: reallocate_im !<
470
471	INTEGER(isp), INTENT(IN) :: m !<
472	INTEGER(isp), INTENT(IN) :: n !<
473
474	INTEGER(isp) :: mold !<
475	INTEGER(isp) :: nold !<
476	INTEGER(isp) :: ierr !<
477
478	ALLOCATE(reallocate_im(n,m),stat=ierr)
479
480	IF (ierr /= 0) CALL &
481	ran_error('reallocate_im: problem in attempt to allocate memory')
482
483	IF (.not. associated(p)) RETURN
484
485	nold = size(p,1)
486	mold = size(p,2)
487
488	reallocate_im(1:min(nold,n),1:min(mold,m)) = &
489	p(1:min(nold,n),1:min(mold,m))
490
491	DEALLOCATE(p)
492
493	END FUNCTION reallocate_im
494
495	!------------------------------------------------------------------------------!
496	! Description:
497	! ------------
498	!> Reallocates the generators state space "ranseeds" to vectors of size length.
499	!------------------------------------------------------------------------------!
500	FUNCTION arth_i(first,increment,n)
501
502	INTEGER(isp), INTENT(IN) :: first !<
503	INTEGER(isp), INTENT(IN) :: increment !<
504	INTEGER(isp), INTENT(IN) :: n !<
505
506	INTEGER(isp), DIMENSION(n) :: arth_i !<
507
508	INTEGER(isp) :: k !<
509	INTEGER(isp) :: k2 !<
510	INTEGER(isp) :: temp !<
511
512	INTEGER(isp), PARAMETER :: npar_arth=16 !<
513	INTEGER(isp), PARAMETER :: npar2_arth=8 !<
514
515	IF (n > 0) arth_i(1) = first
516
517	IF (n <= npar_arth) THEN
518
519	DO k=2,n
520	arth_i(k) = arth_i(k-1)+increment
521	END DO
522
523	ELSE
524
525	DO k=2,npar2_arth
526	arth_i(k) = arth_i(k-1) + increment
527	END DO
528
529	temp = increment * npar2_arth
530	k = npar2_arth
531
532	DO
533	IF (k >= n) EXIT
534	k2 = k + k
535	arth_i(k+1:min(k2,n)) = temp + arth_i(1:min(k,n-k))
536	temp = temp + temp
537	k = k2
538	END DO
539
540	END IF
541
542	END FUNCTION arth_i
543
544	END MODULE random_generator_parallel

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |