Home

Context Navigation

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

Last change on this file since 1695 was 1683, checked in by knoop, 9 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 16.1 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |