Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |