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

Rev	Line
[1850]	1	!> @file random_generator_parallel_mod.f90
[1400]	2	!------------------------------------------------------------------------------!
[2696]	3	! This file is part of the PALM model system.
[1400]	4	!
[2000]	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.
[1400]	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	!
[4360]	17	! Copyright 1997-2020 Leibniz Universitaet Hannover
[1400]	18	!------------------------------------------------------------------------------!
	19	!
	20	! Current revisions:
	21	! -----------------
	22	!
[2173]	23	!
[1400]	24	! Former revisions:
	25	! -----------------
	26	! $Id: random_generator_parallel_mod.f90 4545 2020-05-22 13:17:57Z suehring $
[4545]	27	! Add generator (using parallel mode) returning gaussian distributed random
	28	! number
	29	!
	30	! 4438 2020-03-03 20:49:28Z raasch
[4438]	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
[4182]	37	! Corrected "Former revisions" section
	38	!
	39	! 3655 2019-01-07 16:51:22Z knoop
[3241]	40	! unused variable removed
[4182]	41	!
	42	! 1400 2014-05-09 14:03:54Z knoop
	43	! Initial revision
[3241]	44	!
[1400]	45	!
	46	! Description:
	47	! ------------
[1682]	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).
[4545]	51	!> Moreover, it contains a routine returning a normally distributed random number
	52	!> with mean zero and unity standard deviation.
[1682]	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"
[1400]	62	!------------------------------------------------------------------------------!
[1682]	63	MODULE random_generator_parallel
[1400]	64
	65	USE kinds
[4438]	66
[1400]	67	IMPLICIT NONE
[4438]	68
[1682]	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 !<
[4438]	77
[1682]	78	INTEGER(isp), DIMENSION(:, :), POINTER, SAVE :: ranseeds !<
[4438]	79
[1682]	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 !<
[4438]	86
[1682]	87	INTEGER(isp), DIMENSION(:,:), ALLOCATABLE :: id_random_array !<
	88	INTEGER(isp), DIMENSION(:,:,:), ALLOCATABLE :: seq_random_array !<
[4438]	89
[1682]	90	REAL(wp), SAVE :: amm !<
[4438]	91
[1682]	92	REAL(wp) :: random_dummy=0.0 !<
[4438]	93
[2172]	94	INTERFACE init_parallel_random_generator
[4438]	95	MODULE PROCEDURE init_parallel_random_generator_1d
	96	MODULE PROCEDURE init_parallel_random_generator_2d
[2172]	97	END INTERFACE
[4438]	98
[1400]	99	INTERFACE random_number_parallel
	100	MODULE PROCEDURE ran0_s
	101	END INTERFACE
[4438]	102
[4545]	103	INTERFACE random_number_parallel_gauss
	104	MODULE PROCEDURE gasdev_s
	105	END INTERFACE
	106
[1400]	107	INTERFACE random_seed_parallel
	108	MODULE PROCEDURE random_seed_parallel
	109	END INTERFACE
[4438]	110
[1400]	111	INTERFACE ran_hash
	112	MODULE PROCEDURE ran_hash_v
	113	END INTERFACE
[4438]	114
[1400]	115	INTERFACE reallocate
	116	MODULE PROCEDURE reallocate_iv,reallocate_im
	117	END INTERFACE
[4438]	118
[1400]	119	INTERFACE arth
	120	MODULE PROCEDURE arth_i
	121	END INTERFACE
	122
[4438]	123	PRIVATE
	124
	125	PUBLIC random_number_parallel, random_seed_parallel, random_dummy, &
[4545]	126	id_random_array, seq_random_array, init_parallel_random_generator, &
	127	random_number_parallel_gauss
[4438]	128
[1400]	129	CONTAINS
[4438]	130
[1400]	131	!------------------------------------------------------------------------------!
	132	! Description:
	133	! ------------
[4438]	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	! ------------
[2172]	172	!> Initialize the parallel random number generator for a specific subdomain
	173	!------------------------------------------------------------------------------!
[4438]	174	SUBROUTINE init_parallel_random_generator_2d( nx_l, nys_l, nyn_l, nxl_l, nxr_l )
[2172]	175
	176	USE kinds
	177
	178	USE control_parameters, &
	179	ONLY: ensemble_member_nr
	180
	181	IMPLICIT NONE
	182
[4438]	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
[2172]	188
[4438]	189	INTEGER(iwp) :: i !< grid index x-direction
	190	INTEGER(iwp) :: j !< grid index y-direction
[2172]	191
	192	!-- Allocate ID-array and state-space-array
[4438]	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) )
[2172]	195	seq_random_array = 0
	196	id_random_array = 0
	197
[2696]	198	!-- Asigning an ID to every vertical gridpoint column
	199	!-- dependig on the ensemble run number.
[4438]	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 * &
[2696]	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
[4438]	209	DO i = nxl_l, nxr_l
	210	DO j = nys_l, nyn_l
[2696]	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
[2172]	216
[4438]	217	END SUBROUTINE init_parallel_random_generator_2d
[2172]	218
	219	!------------------------------------------------------------------------------!
	220	! Description:
	221	! ------------
[1682]	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.
[1400]	228	!------------------------------------------------------------------------------!
	229	SUBROUTINE ran0_s(harvest)
	230
	231	USE kinds
	232
	233	IMPLICIT NONE
	234
[1682]	235	REAL(wp), INTENT(OUT) :: harvest !<
[1400]	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	! ------------
[4545]	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	! ------------
[1682]	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.
[1400]	314	!------------------------------------------------------------------------------!
	315	SUBROUTINE ran_init( length )
	316
	317	USE kinds
	318
	319	IMPLICIT NONE
	320
[1682]	321	INTEGER(isp), INTENT(IN) :: length !<
[1400]	322
[2144]	323	INTEGER(isp) :: hg !<
	324	INTEGER(isp) :: hgm !<
	325	INTEGER(isp) :: hgng !<
[1400]	326
[1682]	327	INTEGER(isp) :: new !<
	328	INTEGER(isp) :: j !<
	329	INTEGER(isp) :: hgt !<
[1400]	330
	331	IF ( length < lenran ) RETURN !- Simply return if enough space is already allocated.
	332
[2144]	333	hg=huge(1_isp)
	334	hgm=-hg
	335	hgng=hgm-1
[1400]	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
[2255]	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')
[1400]	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) &
[2255]	367	CALL ran_error('arithmetic assumption 10 failed')
[1400]	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	! ------------
[1682]	413	!> User interface to release the workspace used by the random number routines.
[1400]	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	! ------------
[1682]	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.
[1400]	434	!------------------------------------------------------------------------------!
	435	SUBROUTINE random_seed_parallel( random_sequence, state_size, put, get )
	436
	437	IMPLICIT NONE
	438
[1682]	439	INTEGER(isp), OPTIONAL, INTENT(IN) :: random_sequence !<
	440	INTEGER(isp), OPTIONAL, INTENT(OUT) :: state_size !<
[1400]	441
[1682]	442	INTEGER(isp), DIMENSION(:), OPTIONAL, INTENT(IN) :: put !<
	443	INTEGER(isp), DIMENSION(:), OPTIONAL, INTENT(OUT) :: get !<
[1400]	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	! ------------
[1682]	485	!> DES-like hashing of two 32-bit integers, using shifts,
	486	!> xor's, and adds to make the internal nonlinear function.
[1400]	487	!------------------------------------------------------------------------------!
	488	SUBROUTINE ran_hash_v( il, ir )
	489
	490	IMPLICIT NONE
	491
[1682]	492	INTEGER(isp), DIMENSION(:), INTENT(INOUT) :: il !<
	493	INTEGER(isp), DIMENSION(:), INTENT(INOUT) :: ir !<
[1400]	494
[1682]	495	INTEGER(isp), DIMENSION(size(il)) :: is !<
[1400]	496
[1682]	497	INTEGER(isp) :: j !<
[1400]	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	! ------------
[1682]	514	!> User interface to process error-messages
	515	!> produced by the random_number_parallel module
[1400]	516	!------------------------------------------------------------------------------!
	517	SUBROUTINE ran_error(string)
[2255]	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
[1400]	527	END SUBROUTINE ran_error
	528
	529	!------------------------------------------------------------------------------!
	530	! Description:
	531	! ------------
[1682]	532	!> Reallocates the generators state space "ranseeds" to vectors of size length.
[1400]	533	!------------------------------------------------------------------------------!
	534	FUNCTION reallocate_iv( p, n )
[2255]	535
	536	USE control_parameters, &
	537	ONLY: message_string
[1400]	538
[1682]	539	INTEGER(isp), DIMENSION(:), POINTER :: p !<
	540	INTEGER(isp), DIMENSION(:), POINTER :: reallocate_iv !<
[1400]	541
[1682]	542	INTEGER(isp), INTENT(IN) :: n !<
[1400]	543
[1682]	544	INTEGER(isp) :: nold !<
	545	INTEGER(isp) :: ierr !<
[1400]	546
	547	ALLOCATE(reallocate_iv(n),stat=ierr)
	548
[2255]	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
[1400]	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 )
[2255]	565
	566	USE control_parameters, &
	567	ONLY: message_string
[1400]	568
[1682]	569	INTEGER(isp), DIMENSION(:,:), POINTER :: p !<
	570	INTEGER(isp), DIMENSION(:,:), POINTER :: reallocate_im !<
[1400]	571
[1682]	572	INTEGER(isp), INTENT(IN) :: m !<
	573	INTEGER(isp), INTENT(IN) :: n !<
[1400]	574
[1682]	575	INTEGER(isp) :: mold !<
	576	INTEGER(isp) :: nold !<
	577	INTEGER(isp) :: ierr !<
[1400]	578
	579	ALLOCATE(reallocate_im(n,m),stat=ierr)
	580
[2255]	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
[1400]	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	! ------------
[1682]	601	!> Reallocates the generators state space "ranseeds" to vectors of size length.
[1400]	602	!------------------------------------------------------------------------------!
	603	FUNCTION arth_i(first,increment,n)
	604
[1682]	605	INTEGER(isp), INTENT(IN) :: first !<
	606	INTEGER(isp), INTENT(IN) :: increment !<
	607	INTEGER(isp), INTENT(IN) :: n !<
[1400]	608
[1682]	609	INTEGER(isp), DIMENSION(n) :: arth_i !<
[1400]	610
[1682]	611	INTEGER(isp) :: k !<
	612	INTEGER(isp) :: k2 !<
	613	INTEGER(isp) :: temp !<
[1400]	614
[1682]	615	INTEGER(isp), PARAMETER :: npar_arth=16 !<
	616	INTEGER(isp), PARAMETER :: npar2_arth=8 !<
[1400]	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 |