Home

Context Navigation

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

Last change on this file since 4180 was 4180, checked in by scharf, 5 years ago
removed comments in 'Former revisions' section that are older than 01.01.2019
Property svn:keywords set to `Id`
File size: 18.2 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |