Home

Context Navigation

source: palm/trunk/SOURCE/random_generator_parallel.f90 @ 1504

Last change on this file since 1504 was 1401, checked in by knoop, 11 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 16.0 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |