Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |