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

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

Parallel random number generator added (preliminary version).

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