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

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