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

Last change on this file since 4180 was 4180, checked in by scharf, 22 months ago

removed comments in 'Former revisions' section that are older than 01.01.2019

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