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

Last change on this file since 2145 was 2145, checked in by knoop, 4 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!------------------------------------------------------------------------------!
60MODULE 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  (amm*hgng >= 1.0 .or. amm*hgng <= 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
486END MODULE random_generator_parallel
Note: See TracBrowser for help on using the repository browser.