Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |