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

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