Home

Context Navigation

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

Last change on this file since 2000 was 2000, checked in by knoop, 8 years ago
Forced header and separation lines into 80 columns
Property svn:keywords set to `Id`
File size: 16.2 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |