Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |