Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |