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

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

Implemented PALM specific error message handling

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