Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |