Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |