Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |