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

Last change on this file since 1401 was 1401, checked in by knoop, 7 years ago

last commit documented

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