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

Last change on this file since 2173 was 2173, checked in by knoop, 4 years ago

last commit documented

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