Home

Context Navigation

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

Last change on this file since 2223 was 2173, checked in by knoop, 8 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 18.4 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |