Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |