Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |