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

Last change on this file since 2703 was 2696, checked in by kanani, 7 years ago

Merge of branch palm4u into trunk

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