Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |