source: palm/trunk/SOURCE/random_generator_parallel.f90 @ 1682

Last change on this file since 1682 was 1682, checked in by knoop, 9 years ago

Code annotations made doxygen readable

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