Home

Context Navigation

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

Last change on this file since 1682 was 1682, checked in by knoop, 9 years ago
Code annotations made doxygen readable
Property svn:keywords set to `Id`
File size: 16.1 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |