Home

Context Navigation

random_generator_parallel_mod.f90 @ 4892

Last change on this file since 4892 was 4828, checked in by Giersch, 4 years ago
Copyright updated to year 2021, interface pmc_sort removed to accelarate the nesting code
Property svn:keywords set to `Id`
File size: 22.4 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 terms of the GNU General
6	! Public License as published by the Free Software Foundation, either version 3 of the License, or
7	! (at your option) any later version.
8	!
9	! PALM is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the
10	! implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
11	! Public License for more details.
12	!
13	! You should have received a copy of the GNU General Public License along with PALM. If not, see
14	! <http://www.gnu.org/licenses/>.
15	!
16	! Copyright 1997-2021 Leibniz Universitaet Hannover
17	!--------------------------------------------------------------------------------------------------!
18	!
19	!
20	! Current revisions:
21	! -----------------
22	!
23	!
24	! Former revisions:
25	! -----------------
26	! $Id: random_generator_parallel_mod.f90 4828 2021-01-05 11:21:41Z suehring $
27	! File re-formatted to follow the PALM coding standard
28	!
29	! 4545 2020-05-22 13:17:57Z schwenkel
30	! Add generator (using parallel mode) returning gaussian distributed random number
31	!
32	! 4438 2020-03-03 20:49:28Z raasch
33	! - Rename variables to avoid confusion with subdomain grid indices
34	! - Some formatting adjustments
35	! - New routine to initialize spatial 1D-arrays independent on the 2D random-number array
36	!
37	! 4360 2020-01-07 11:25:50Z suehring
38	! Corrected "Former revisions" section
39	!
40	! 3655 2019-01-07 16:51:22Z knoop
41	! Unused variable removed
42	!
43	! 1400 2014-05-09 14:03:54Z knoop
44	! Initial revision
45	!
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 (exclusive of the end point
52	!> values). Moreover, it contains a routine returning a normally distributed random number with mean
53	!> zero and unity standard deviation. Additionally, it provides the generator with five integers for
54	!> use as initial state space. The first tree integers (iran, jran, kran) are maintained as non
55	!> negative values, while the last two (mran, nran) have 32-bit nonzero values. Also provided by
56	!> this module is support for initializing or reinitializing the state space to a desired standard
57	!> sequence number, hashing the initial values to random values and allocating and deallocating the
58	!> internal workspace.
59	!>
60	!> This routine is taken from the "numerical recipies vol. 2"
61	!--------------------------------------------------------------------------------------------------!
62	MODULE random_generator_parallel
63
64	USE kinds
65
66	IMPLICIT NONE
67
68	INTEGER(isp), SAVE :: lenran = 0 !<
69	INTEGER(isp), SAVE :: seq = 0 !<
70	INTEGER(isp), SAVE :: iran0 !<
71	INTEGER(isp), SAVE :: jran0 !<
72	INTEGER(isp), SAVE :: kran0 !<
73	INTEGER(isp), SAVE :: mran0 !<
74	INTEGER(isp), SAVE :: nran0 !<
75	INTEGER(isp), SAVE :: rans !<
76
77	INTEGER(isp), DIMENSION(:, :), POINTER, SAVE :: ranseeds !<
78
79	INTEGER(isp), DIMENSION(:), POINTER, SAVE :: iran !<
80	INTEGER(isp), DIMENSION(:), POINTER, SAVE :: jran !<
81	INTEGER(isp), DIMENSION(:), POINTER, SAVE :: kran !<
82	INTEGER(isp), DIMENSION(:), POINTER, SAVE :: mran !<
83	INTEGER(isp), DIMENSION(:), POINTER, SAVE :: nran !<
84	INTEGER(isp), DIMENSION(:), POINTER, SAVE :: ranv !<
85
86	INTEGER(isp), DIMENSION(:,:), ALLOCATABLE :: id_random_array !<
87	INTEGER(isp), DIMENSION(:,:,:), ALLOCATABLE :: seq_random_array !<
88
89	REAL(wp), SAVE :: amm !<
90
91	REAL(wp) :: random_dummy = 0.0 !<
92
93	INTERFACE init_parallel_random_generator
94	MODULE PROCEDURE init_parallel_random_generator_1d
95	MODULE PROCEDURE init_parallel_random_generator_2d
96	END INTERFACE
97
98	INTERFACE random_number_parallel
99	MODULE PROCEDURE ran0_s
100	END INTERFACE
101
102	INTERFACE random_number_parallel_gauss
103	MODULE PROCEDURE gasdev_s
104	END INTERFACE
105
106	INTERFACE random_seed_parallel
107	MODULE PROCEDURE random_seed_parallel
108	END INTERFACE
109
110	INTERFACE ran_hash
111	MODULE PROCEDURE ran_hash_v
112	END INTERFACE
113
114	INTERFACE reallocate
115	MODULE PROCEDURE reallocate_iv,reallocate_im
116	END INTERFACE
117
118	INTERFACE arth
119	MODULE PROCEDURE arth_i
120	END INTERFACE
121
122	PRIVATE
123
124	PUBLIC random_number_parallel, &
125	random_seed_parallel, &
126	random_dummy, &
127	id_random_array, &
128	seq_random_array, &
129	init_parallel_random_generator, &
130	random_number_parallel_gauss
131
132	CONTAINS
133
134	!--------------------------------------------------------------------------------------------------!
135	! Description:
136	! ------------
137	!> Initialize the parallel random number generator for a 1-dimensional array.
138	!--------------------------------------------------------------------------------------------------!
139	SUBROUTINE init_parallel_random_generator_1d( nxy, ns, ne, id_rand, seq_rand )
140
141	USE control_parameters, &
142	ONLY: ensemble_member_nr
143
144	INTEGER(isp), INTENT(IN) :: nxy !< constant scaling with grid dimensions
145	INTEGER(isp), INTENT(IN) :: ns !< start index on subdomain
146	INTEGER(isp), INTENT(IN) :: ne !< end index on subdomain
147
148	INTEGER(iwp) :: j !< loop index
149
150	INTEGER(isp), DIMENSION(ns:ne) :: id_rand !< initial IDs
151	INTEGER(isp), DIMENSION(5,ns:ne) :: seq_rand !< initial random seeds
152
153	!
154	!-- Asigning an ID to every vertical gridpoint column dependig on the ensemble run number.
155	DO j = ns, ne
156	id_rand(j) = j * ( nxy + 1.0_wp ) + 1.0_wp + 1E6 * ( ensemble_member_nr - 1000 )
157	ENDDO
158	!
159	!-- Initializing with random_seed_parallel for every vertical gridpoint column.
160	random_dummy = 0
161	DO j = ns, ne
162	CALL random_seed_parallel( random_sequence=id_rand(j) )
163	CALL random_number_parallel( random_dummy )
164	CALL random_seed_parallel( get=seq_rand(:,j) )
165	ENDDO
166
167	END SUBROUTINE init_parallel_random_generator_1d
168
169	!--------------------------------------------------------------------------------------------------!
170	! Description:
171	! ------------
172	!> Initialize the parallel random number generator for a specific subdomain.
173	!--------------------------------------------------------------------------------------------------!
174	SUBROUTINE init_parallel_random_generator_2d( nx_l, nys_l, nyn_l, nxl_l, nxr_l )
175
176	USE kinds
177
178	USE control_parameters, &
179	ONLY: ensemble_member_nr
180
181	IMPLICIT NONE
182
183	INTEGER(isp), INTENT(IN) :: nx_l !< constant
184	INTEGER(isp), INTENT(IN) :: nys_l !< local lower subdomain bound index in y-direction
185	INTEGER(isp), INTENT(IN) :: nyn_l !< local upper subdomain bound index in y-direction
186	INTEGER(isp), INTENT(IN) :: nxl_l !< local lower subdomain bound index in x-direction
187	INTEGER(isp), INTENT(IN) :: nxr_l !< local upper subdomain bound index in x-direction
188
189	INTEGER(iwp) :: i !< grid index x-direction
190	INTEGER(iwp) :: j !< grid index y-direction
191
192	!-- Allocate ID-array and state-space-array
193	ALLOCATE ( seq_random_array(5,nys_l:nyn_l,nxl_l:nxr_l) )
194	ALLOCATE ( id_random_array(nys_l:nyn_l,nxl_l:nxr_l) )
195	seq_random_array = 0
196	id_random_array = 0
197
198	!-- Asigning an ID to every vertical gridpoint column dependig on the ensemble run number.
199	DO i = nxl_l, nxr_l
200	DO j = nys_l, nyn_l
201	id_random_array(j,i) = j * ( nx_l + 1.0_wp ) + i + 1.0_wp + 1E6 * &
202	( ensemble_member_nr - 1000 )
203	ENDDO
204	ENDDO
205	!-- Initializing with random_seed_parallel for every vertical gridpoint column.
206	random_dummy = 0
207	DO i = nxl_l, nxr_l
208	DO j = nys_l, nyn_l
209	CALL random_seed_parallel( random_sequence = id_random_array(j,i) )
210	CALL random_number_parallel( random_dummy )
211	CALL random_seed_parallel( get = seq_random_array(:,j,i) )
212	ENDDO
213	ENDDO
214
215	END SUBROUTINE init_parallel_random_generator_2d
216
217	!--------------------------------------------------------------------------------------------------!
218	! Description:
219	! ------------
220	!> Lagged Fibonacci generator combined with a Marsaglia shift sequence. Returns as harvest a uniform
221	!> random deviate between 0.0 and 1.0 (exclusive of the end point values). This generator has the
222	!> same calling and initialization conventions as Fortran 90's random_number routine.
223	!> Use random_seed_parallel to initialize or reinitialize to a particular sequence.
224	!> The period of this generator is about 2.0 x 10^28, and it fully vectorizes.
225	!> Validity of the integer model assumed by this generator is tested at initialization.
226	!--------------------------------------------------------------------------------------------------!
227	SUBROUTINE ran0_s( harvest )
228
229	USE kinds
230
231	IMPLICIT NONE
232
233	REAL(wp), INTENT(OUT) :: harvest !<
234	!
235	!-- Initialization routine in ran_state.
236	IF ( lenran < 1 ) CALL ran_init( 1 )
237	!
238	!-- Update Fibonacci generator, which has period p^2 + p + 1, p = 2^31 - 69.
239	rans = iran0 - kran0
240
241	IF ( rans < 0 ) rans = rans + 2147483579_isp
242
243	iran0 = jran0
244	jran0 = kran0
245	kran0 = rans
246	!
247	!-- Update Marsaglia shift sequence with period 2^32 - 1.
248	nran0 = IEOR( nran0, ISHFT( nran0, 13 ) )
249	nran0 = IEOR( nran0, ISHFT( nran0, -17 ) )
250	nran0 = IEOR( nran0, ISHFT( nran0, 5 ) )
251	!
252	!-- Combine the generators.
253	rans = IEOR( nran0, rans )
254	!
255	!-- Make the result positive definite (note that amm is negative).
256	harvest = amm * MERGE( rans, NOT( rans ), rans < 0 )
257
258	END SUBROUTINE ran0_s
259
260	!--------------------------------------------------------------------------------------------------!
261	! Description:
262	! ------------
263	!> Returns in harvest a normally distributed deviate with zero mean and unit variance, using ran0_s
264	!> as the source of uniform deviates. Following Numerical Recipes in Fortran90 (Press et al., 2nd
265	!> edition, 1996, pp 1152ff). Note, instead of ran1_s, ran0_s is used.
266	!--------------------------------------------------------------------------------------------------!
267	SUBROUTINE gasdev_s( harvest )
268
269	REAL(wp), INTENT(OUT) :: harvest !<
270
271	REAL(wp) :: rsq !<
272	REAL(wp) :: v1 !<
273	REAL(wp) :: v2 !<
274	REAL(wp), SAVE :: g !<
275
276	LOGICAL, SAVE :: gaus_stored = .FALSE. !<
277	!
278	!-- We have an extra deviate handy, so return it, and unset the flag.
279	IF ( gaus_stored ) THEN
280	harvest = g
281	gaus_stored = .FALSE.
282	!
283	!-- We donât have an extra deviate handy, so pick two uniform numbers in the square extending
284	!-- from -1 to +1 in each direction
285	ELSE
286	DO
287	CALL ran0_s( v1 )
288	CALL ran0_s( v2 )
289	v1 = 2.0_wp * v1 - 1.0_wp
290	v2 = 2.0_wp * v2 - 1.0_wp
291	!
292	!-- See if they are in the unit circle
293	rsq = v12 + v22
294	!
295	!-- Otherwise try again.
296	IF ( rsq > 0.0 .AND. rsq < 1.0 ) EXIT
297	ENDDO
298	!
299	!-- Now make the Box-Muller transformation to get two normal deviates.
300	!-- Return one and save the other for next time. Set flag.
301	rsq = SQRT( - 2.0_sp * LOG( rsq ) / rsq )
302	harvest = v1 * rsq
303	g = v2 * rsq
304	gaus_stored = .TRUE.
305	ENDIF
306
307	END SUBROUTINE gasdev_s
308
309	!--------------------------------------------------------------------------------------------------!
310	! Description:
311	! ------------
312	!> Initialize or reinitialize the random generator state space to vectors of size length.
313	!> The saved variable seq is hashed (via calls to the module routine ran_hash) to create unique
314	!> starting seeds, different for each vector component.
315	!--------------------------------------------------------------------------------------------------!
316	SUBROUTINE ran_init( length )
317
318	USE kinds
319
320	IMPLICIT NONE
321
322	INTEGER(isp), INTENT(IN) :: length !<
323
324	INTEGER(isp) :: hg !<
325	INTEGER(isp) :: hgm !<
326	INTEGER(isp) :: hgng !<
327
328	INTEGER(isp) :: new !<
329	INTEGER(isp) :: j !<
330	INTEGER(isp) :: hgt !<
331	!
332	!-- Simply return if enough space is already allocated.
333	IF ( length < lenran ) RETURN
334
335	hg = HUGE( 1_isp )
336	hgm = - hg
337	hgng = hgm - 1
338	hgt = hg
339
340	!-- The following lines check that kind value isp is in fact a 32-bit integer with the usual
341	!-- properties that we expect it to have (under negation and wrap-around addition).
342	!-- If all of these tests are satisfied, then the routines that use this module are portable, even
343	!-- though they go beyond Fortran 90's integer model.
344
345	IF ( hg /= 2147483647 ) CALL ran_error( 'arithmetic assumption 1 failed' )
346	IF ( hgng >= 0 ) CALL ran_error( 'arithmetic assumption 2 failed' )
347	IF ( hgt + 1 /= hgng ) CALL ran_error( 'arithmetic assumption 3 failed' )
348	IF ( NOT( hg ) >= 0 ) CALL ran_error( 'arithmetic assumption 4 failed' )
349	IF ( NOT( hgng ) < 0 ) CALL ran_error( 'arithmetic assumption 5 failed' )
350	IF ( hg + hgng >= 0 ) CALL ran_error( 'arithmetic assumption 6 failed' )
351	IF ( NOT( - 1_isp ) < 0 ) CALL ran_error( 'arithmetic assumption 7 failed' )
352	IF ( NOT( 0_isp ) >= 0 ) CALL ran_error( 'arithmetic assumption 8 failed' )
353	IF ( NOT( 1_isp ) >= 0 ) CALL ran_error( 'arithmetic assumption 9 failed' )
354
355	IF ( lenran > 0 ) THEN ! Reallocate space, or ...
356
357	ranseeds => reallocate( ranseeds, length, 5 )
358	ranv => reallocate( ranv, length - 1 )
359	new = lenran + 1
360
361	ELSE ! Allocate space.
362
363	ALLOCATE( ranseeds(length,5) )
364	ALLOCATE( ranv(length-1) )
365	new = 1 !- Index of first location not yet initialized.
366	amm = NEAREST( 1.0_wp, - 1.0_wp ) / hgng
367	!
368	!-- Use of nearest is to ensure that returned random deviates are strictly less than 1.0.
369	IF ( amm * hgng >= 1.0 .OR. amm * hgng <= 0.0 ) &
370	CALL ran_error( 'arithmetic assumption 10 failed' )
371
372	END IF
373	!
374	!-- Set starting values, unique by seq and vector component.
375	ranseeds(new:,1) = seq
376	ranseeds(new:,2:5) = SPREAD( arth( new, 1, SIZE( ranseeds(new:,1) ) ), 2, 4 )
377
378	DO j = 1, 4 ! Hash them.
379	CALL ran_hash( ranseeds(new:,j), ranseeds(new:,j+1) )
380	END DO
381
382	WHERE ( ranseeds(new:,1:3) < 0 ) &
383	ranseeds(new:,1:3) = NOT( ranseeds(new:,1:3) ) ! Enforce nonnegativity.
384
385	WHERE ( ranseeds(new:,4:5) == 0 ) ranseeds(new:,4:5) = 1 ! Enforce nonzero.
386
387	IF ( new == 1 ) THEN ! Set scalar seeds.
388
389	iran0 = ranseeds(1,1)
390	jran0 = ranseeds(1,2)
391	kran0 = ranseeds(1,3)
392	mran0 = ranseeds(1,4)
393	nran0 = ranseeds(1,5)
394	rans = nran0
395
396	END IF
397
398	IF ( length > 1 ) THEN ! Point to vector seeds.
399
400	iran => ranseeds(2:,1)
401	jran => ranseeds(2:,2)
402	kran => ranseeds(2:,3)
403	mran => ranseeds(2:,4)
404	nran => ranseeds(2:,5)
405	ranv = nran
406
407	END IF
408
409	lenran = length
410
411	END SUBROUTINE ran_init
412
413	!--------------------------------------------------------------------------------------------------!
414	! Description:
415	! ------------
416	!> User interface to release the workspace used by the random number routines.
417	!--------------------------------------------------------------------------------------------------!
418	SUBROUTINE ran_deallocate
419
420	IF ( lenran > 0 ) THEN
421
422	DEALLOCATE( ranseeds, ranv )
423	NULLIFY( ranseeds, ranv, iran, jran, kran, mran, nran )
424	lenran = 0
425
426	END IF
427
428	END SUBROUTINE ran_deallocate
429
430	!--------------------------------------------------------------------------------------------------!
431	! Description:
432	! ------------
433	!> User interface for seeding the random number routines.
434	!> Syntax is exactly like Fortran 90's random_seed routine, with one additional argument keyword:
435	!> random_sequence, set to any integer value, causes an immediate new initialization, seeded by that
436	!> integer.
437	!--------------------------------------------------------------------------------------------------!
438	SUBROUTINE random_seed_parallel( random_sequence, state_size, put, get )
439
440	IMPLICIT NONE
441
442	INTEGER(isp), OPTIONAL, INTENT(IN) :: random_sequence !<
443	INTEGER(isp), OPTIONAL, INTENT(OUT) :: state_size !<
444
445	INTEGER(isp), DIMENSION(:), OPTIONAL, INTENT(IN) :: put !<
446	INTEGER(isp), DIMENSION(:), OPTIONAL, INTENT(OUT) :: get !<
447
448	IF ( PRESENT( state_size ) ) THEN
449
450	state_size = 5 * lenran
451
452	ELSE IF ( PRESENT( put ) ) THEN
453
454	IF ( lenran == 0 ) RETURN
455
456	ranseeds = RESHAPE( put, SHAPE( ranseeds ) )
457	!
458	!-- Enforce nonnegativity and nonzero conditions on any user-supplied seeds.
459	WHERE ( ranseeds(:,1:3) < 0 ) ranseeds(:,1:3) = NOT( ranseeds(:,1:3) )
460
461	WHERE ( ranseeds(:,4:5) == 0 ) ranseeds(:,4:5) = 1
462
463	iran0 = ranseeds(1,1)
464	jran0 = ranseeds(1,2)
465	kran0 = ranseeds(1,3)
466	mran0 = ranseeds(1,4)
467	nran0 = ranseeds(1,5)
468
469	ELSE IF ( PRESENT( get ) ) THEN
470
471	IF ( lenran == 0 ) RETURN
472
473	ranseeds(1,1:5) = (/ iran0, jran0, kran0, mran0, nran0 /)
474	get = RESHAPE( ranseeds, SHAPE( get ) )
475
476	ELSE IF ( PRESENT( random_sequence ) ) THEN
477
478	CALL ran_deallocate
479	seq = random_sequence
480
481	END IF
482
483	END SUBROUTINE random_seed_parallel
484
485	!--------------------------------------------------------------------------------------------------!
486	! Description:
487	! ------------
488	!> DES-like hashing of two 32-bit integers, using shifts, xor's, and adds to make the internal
489	!> nonlinear function.
490	!--------------------------------------------------------------------------------------------------!
491	SUBROUTINE ran_hash_v( il, ir )
492
493	IMPLICIT NONE
494
495	INTEGER(isp), DIMENSION(:), INTENT(INOUT) :: il !<
496	INTEGER(isp), DIMENSION(:), INTENT(INOUT) :: ir !<
497
498	INTEGER(isp), DIMENSION(size(il)) :: is !<
499
500	INTEGER(isp) :: j !<
501
502	DO j = 1, 4
503
504	is = ir
505	ir = IEOR( ir, ISHFT( ir, 5 ) ) + 1422217823
506	ir = IEOR( ir, ISHFT( ir, - 16 ) ) + 1842055030
507	ir = IEOR( ir, ISHFT( ir, 9 ) ) + 80567781
508	ir = IEOR( il, ir )
509	il = is
510	END DO
511
512	END SUBROUTINE ran_hash_v
513
514	!--------------------------------------------------------------------------------------------------!
515	! Description:
516	! ------------
517	!> User interface to process error-messages produced by the random_number_parallel module.
518	!--------------------------------------------------------------------------------------------------!
519	SUBROUTINE ran_error( string )
520
521	USE control_parameters, &
522	ONLY: message_string
523
524	CHARACTER(LEN=*), INTENT(IN) :: string !< Error message string
525
526	message_string = 'incompatible integer arithmetic: ' // TRIM( string )
527	CALL message( 'ran_init', 'PA0453', 1, 2, 0, 6, 0 )
528
529	END SUBROUTINE ran_error
530
531	!--------------------------------------------------------------------------------------------------!
532	! Description:
533	! ------------
534	!> Reallocates the generators state space "ranseeds" to vectors of size length.
535	!--------------------------------------------------------------------------------------------------!
536	FUNCTION reallocate_iv( p, n )
537
538	USE control_parameters, &
539	ONLY: message_string
540
541	INTEGER(isp), DIMENSION(:), POINTER :: p !<
542	INTEGER(isp), DIMENSION(:), POINTER :: reallocate_iv !<
543
544	INTEGER(isp), INTENT(IN) :: n !<
545
546	INTEGER(isp) :: nold !<
547	INTEGER(isp) :: ierr !<
548
549	ALLOCATE( reallocate_iv(n), stat = ierr )
550
551	IF ( ierr /= 0 ) THEN
552	message_string = 'problem in attempt to allocate memory'
553	CALL message( 'reallocate_iv', 'PA0454', 1, 2, 0, 6, 0 )
554	END IF
555
556	IF ( .NOT. ASSOCIATED( p ) ) RETURN
557
558	nold = SIZE( p )
559
560	reallocate_iv(1:MIN( nold, n )) = p(1:MIN( nold, n ))
561
562	DEALLOCATE( p )
563
564	END FUNCTION reallocate_iv
565
566	FUNCTION reallocate_im( p, n, m )
567
568	USE control_parameters, &
569	ONLY: message_string
570
571	INTEGER(isp), DIMENSION(:,:), POINTER :: p !<
572	INTEGER(isp), DIMENSION(:,:), POINTER :: reallocate_im !<
573
574	INTEGER(isp), INTENT(IN) :: m !<
575	INTEGER(isp), INTENT(IN) :: n !<
576
577	INTEGER(isp) :: mold !<
578	INTEGER(isp) :: nold !<
579	INTEGER(isp) :: ierr !<
580
581	ALLOCATE( reallocate_im(n,m), stat = ierr )
582
583	IF ( ierr /= 0 ) THEN
584	message_string = 'problem in attempt to allocate memory'
585	CALL message( 'reallocate_im', 'PA0454', 1, 2, 0, 6, 0 )
586	END IF
587
588	IF ( .NOT. ASSOCIATED( p ) ) RETURN
589
590	nold = SIZE( p, 1 )
591	mold = SIZE( p, 2 )
592
593	reallocate_im(1:MIN( nold, n ),1:MIN( mold, m )) = p(1:MIN( nold, n ),1:MIN( mold, m ))
594
595	DEALLOCATE(p)
596
597	END FUNCTION reallocate_im
598
599	!--------------------------------------------------------------------------------------------------!
600	! Description:
601	! ------------
602	!> Reallocates the generators state space "ranseeds" to vectors of size length.
603	!--------------------------------------------------------------------------------------------------!
604	FUNCTION arth_i( first, increment, n )
605
606	INTEGER(isp), INTENT(IN) :: first !<
607	INTEGER(isp), INTENT(IN) :: increment !<
608	INTEGER(isp), INTENT(IN) :: n !<
609
610	INTEGER(isp), DIMENSION(n) :: arth_i !<
611
612	INTEGER(isp) :: k !<
613	INTEGER(isp) :: k2 !<
614	INTEGER(isp) :: temp !<
615
616	INTEGER(isp), PARAMETER :: npar_arth = 16 !<
617	INTEGER(isp), PARAMETER :: npar2_arth = 8 !<
618
619	IF ( n > 0 ) arth_i(1) = first
620
621	IF ( n <= npar_arth ) THEN
622
623	DO k = 2, n
624	arth_i(k) = arth_i(k-1) + increment
625	END DO
626
627	ELSE
628
629	DO k = 2, npar2_arth
630	arth_i(k) = arth_i(k-1) + increment
631	END DO
632
633	temp = increment * npar2_arth
634	k = npar2_arth
635
636	DO
637	IF ( k >= n ) EXIT
638	k2 = k + k
639	arth_i(k+1:MIN( k2, n )) = temp + arth_i(1:MIN( k, n - k ))
640	temp = temp + temp
641	k = k2
642	END DO
643
644	END IF
645
646	END FUNCTION arth_i
647
648	END MODULE random_generator_parallel

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

Context Navigation

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

Download in other formats: