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

Last change on this file since 4545 was 4545, checked in by schwenkel, 19 months ago

Add gaussian random number generator to parallel random generator and using parallel random number in lpm

  • Property svn:keywords set to Id
File size: 22.4 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-2020 Leibniz Universitaet Hannover
18!------------------------------------------------------------------------------!
19!
20! Current revisions:
21! -----------------
22!
23!
24! Former revisions:
25! -----------------
26! $Id: random_generator_parallel_mod.f90 4545 2020-05-22 13:17:57Z schwenkel $
27! Add generator (using parallel mode) returning gaussian distributed random
28! number
29!
30! 4438 2020-03-03 20:49:28Z raasch
31! - Rename variables to avoid confusion with subdomain grid indices
32! - Some formatting adjustments
33! - New routine to initialize spatial 1D-arrays independent on the 2D random-
34!   number array
35!
36! 4360 2020-01-07 11:25:50Z suehring
37! Corrected "Former revisions" section
38!
39! 3655 2019-01-07 16:51:22Z knoop
40! unused variable removed
41!
42! 1400 2014-05-09 14:03:54Z knoop
43! Initial revision
44!
45!
46! Description:
47! ------------
48!> This module contains and supports the random number generating routine ran_parallel.
49!> ran_parallel returns a uniform random deviate between 0.0 and 1.0
50!> (exclusive of the end point values).
51!> Moreover, it contains a routine returning a normally distributed random number
52!> with mean zero and unity standard deviation.
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   USE kinds
66
67   IMPLICIT NONE
68
69   INTEGER(isp), SAVE :: lenran=0             !<
70   INTEGER(isp), SAVE :: seq=0                !<
71   INTEGER(isp), SAVE :: iran0                !<
72   INTEGER(isp), SAVE :: jran0                !<
73   INTEGER(isp), SAVE :: kran0                !<
74   INTEGER(isp), SAVE :: mran0                !<
75   INTEGER(isp), SAVE :: nran0                !<
76   INTEGER(isp), SAVE :: rans                 !<
77
78   INTEGER(isp), DIMENSION(:, :), POINTER, SAVE :: ranseeds   !<
79
80   INTEGER(isp), DIMENSION(:), POINTER, SAVE :: iran   !<
81   INTEGER(isp), DIMENSION(:), POINTER, SAVE :: jran   !<
82   INTEGER(isp), DIMENSION(:), POINTER, SAVE :: kran   !<
83   INTEGER(isp), DIMENSION(:), POINTER, SAVE :: mran   !<
84   INTEGER(isp), DIMENSION(:), POINTER, SAVE :: nran   !<
85   INTEGER(isp), DIMENSION(:), POINTER, SAVE :: ranv   !<
86
87   INTEGER(isp), DIMENSION(:,:), ALLOCATABLE   ::  id_random_array    !<
88   INTEGER(isp), DIMENSION(:,:,:), ALLOCATABLE ::  seq_random_array   !<
89
90   REAL(wp), SAVE :: amm   !<
91
92   REAL(wp) :: random_dummy=0.0   !<
93
94   INTERFACE init_parallel_random_generator
95      MODULE PROCEDURE init_parallel_random_generator_1d
96      MODULE PROCEDURE init_parallel_random_generator_2d
97   END INTERFACE
98
99   INTERFACE random_number_parallel
100      MODULE PROCEDURE ran0_s
101   END INTERFACE
102
103   INTERFACE random_number_parallel_gauss
104      MODULE PROCEDURE gasdev_s
105   END INTERFACE
106
107   INTERFACE random_seed_parallel
108      MODULE PROCEDURE random_seed_parallel
109   END INTERFACE
110
111   INTERFACE ran_hash
112      MODULE PROCEDURE ran_hash_v
113   END INTERFACE
114
115   INTERFACE reallocate
116      MODULE PROCEDURE reallocate_iv,reallocate_im
117   END INTERFACE
118
119   INTERFACE arth
120      MODULE PROCEDURE arth_i
121   END INTERFACE
122
123   PRIVATE
124
125   PUBLIC random_number_parallel, random_seed_parallel, random_dummy,          &
126          id_random_array, seq_random_array, init_parallel_random_generator,   &
127          random_number_parallel_gauss
128
129 CONTAINS
130
131!------------------------------------------------------------------------------!
132! Description:
133! ------------
134!> Initialize the parallel random number generator for a 1-dimensional array.
135!------------------------------------------------------------------------------!
136   SUBROUTINE init_parallel_random_generator_1d( nxy, ns, ne, id_rand, seq_rand )
137
138      USE control_parameters,                                                  &
139          ONLY: ensemble_member_nr
140
141      INTEGER(isp), INTENT(IN) ::  nxy  !< constant scaling with grid dimensions
142      INTEGER(isp), INTENT(IN) ::  ns   !< start index on subdomain
143      INTEGER(isp), INTENT(IN) ::  ne   !< end index on subdomain
144
145      INTEGER(iwp) ::  j   !< loop index
146
147      INTEGER(isp), DIMENSION(ns:ne)   ::  id_rand  !< initial IDs
148      INTEGER(isp), DIMENSION(5,ns:ne) ::  seq_rand !< initial random seeds
149
150!
151!--   Asigning an ID to every vertical gridpoint column
152!--   dependig on the ensemble run number.
153      DO  j = ns, ne
154         id_rand(j) = j * ( nxy + 1.0_wp ) + 1.0_wp + 1E6 *                    &
155                                ( ensemble_member_nr - 1000 )
156      ENDDO
157!
158!--   Initializing with random_seed_parallel for every vertical
159!--   gridpoint column.
160      random_dummy=0
161      DO  j = ns, ne
162         CALL random_seed_parallel( random_sequence=id_rand(j) )
163         CALL random_number_parallel( random_dummy )
164         CALL random_seed_parallel( get=seq_rand(:,j) )
165      ENDDO
166
167   END SUBROUTINE init_parallel_random_generator_1d
168
169!------------------------------------------------------------------------------!
170! Description:
171! ------------
172!> Initialize the parallel random number generator for a specific subdomain
173!------------------------------------------------------------------------------!
174   SUBROUTINE init_parallel_random_generator_2d( nx_l, nys_l, nyn_l, nxl_l, nxr_l )
175
176      USE kinds
177
178      USE control_parameters,                                                  &
179          ONLY: ensemble_member_nr
180
181      IMPLICIT NONE
182
183      INTEGER(isp), INTENT(IN) ::  nx_l    !< constant
184      INTEGER(isp), INTENT(IN) ::  nys_l   !< local lower subdomain bound index in y-direction
185      INTEGER(isp), INTENT(IN) ::  nyn_l   !< local upper subdomain bound index in y-direction
186      INTEGER(isp), INTENT(IN) ::  nxl_l   !< local lower subdomain bound index in x-direction
187      INTEGER(isp), INTENT(IN) ::  nxr_l   !< local upper subdomain bound index in x-direction
188
189      INTEGER(iwp) ::  i   !< grid index x-direction
190      INTEGER(iwp) ::  j   !< grid index y-direction
191
192!--   Allocate ID-array and state-space-array
193      ALLOCATE ( seq_random_array(5,nys_l:nyn_l,nxl_l:nxr_l) )
194      ALLOCATE ( id_random_array(nys_l:nyn_l,nxl_l:nxr_l) )
195      seq_random_array = 0
196      id_random_array  = 0
197
198!--   Asigning an ID to every vertical gridpoint column
199!--   dependig on the ensemble run number.
200      DO  i = nxl_l, nxr_l
201         DO  j = nys_l, nyn_l
202            id_random_array(j,i) = j * ( nx_l + 1.0_wp ) + i + 1.0_wp + 1E6 *  &
203                                   ( ensemble_member_nr - 1000 )
204         ENDDO
205      ENDDO
206!--   Initializing with random_seed_parallel for every vertical
207!--   gridpoint column.
208      random_dummy=0
209      DO  i = nxl_l, nxr_l
210         DO  j = nys_l, nyn_l
211            CALL random_seed_parallel (random_sequence=id_random_array(j, i))
212            CALL random_number_parallel (random_dummy)
213            CALL random_seed_parallel (get=seq_random_array(:, j, i))
214         ENDDO
215      ENDDO
216
217   END SUBROUTINE init_parallel_random_generator_2d
218 
219!------------------------------------------------------------------------------!
220! Description:
221! ------------
222!> Lagged Fibonacci generator combined with a Marsaglia shiftsequence.
223!> Returns as harvest a uniform random deviate between 0.0 and 1.0 (exclusive of the end point values).
224!> This generator has the same calling and initialization conventions as Fortran 90's random_number routine.
225!> Use random_seed_parallel to initialize or reinitialize to a particular sequence.
226!> The period of this generator is about 2.0 x 10^28, and it fully vectorizes.
227!> Validity of the integer model assumed by this generator is tested at initialization.
228!------------------------------------------------------------------------------!
229   SUBROUTINE ran0_s(harvest)
230
231      USE kinds
232     
233      IMPLICIT NONE
234     
235      REAL(wp), INTENT(OUT) :: harvest   !<
236     
237      IF  (lenran < 1) CALL ran_init(1)  !- Initialization routine in ran_state.
238     
239      !- Update Fibonacci generator, which has period p^2 + p + 1, p = 2^31 - 69.
240      rans = iran0 - kran0   
241     
242      IF  (rans < 0) rans = rans + 2147483579_isp
243     
244      iran0 = jran0
245      jran0 = kran0
246      kran0 = rans
247     
248      nran0 = ieor( nran0, ishft (nran0, 13) ) !- Update Marsaglia shift sequence with period 2^32 - 1.
249      nran0 = ieor( nran0, ishft (nran0, -17) )
250      nran0 = ieor( nran0, ishft (nran0, 5) )
251     
252      rans  = ieor( nran0, rans )   !- Combine the generators.
253     
254      harvest = amm * merge( rans, not(rans), rans < 0 ) !- Make the result positive definite (note that amm is negative).
255     
256   END SUBROUTINE ran0_s
257
258!------------------------------------------------------------------------------!
259! Description:
260! ------------
261!> Returns in harvest a normally distributed deviate with zero mean and unit
262!> variance, using ran0_s as the source of uniform deviates. Following
263!> Numerical Recipes in Fortran90 (Press et al., 2nd edition, 1996, pp 1152ff).
264!> Note, instead of ran1_s, ran0_s is used.
265!------------------------------------------------------------------------------!
266   SUBROUTINE gasdev_s(harvest)
267
268      REAL(wp), INTENT(OUT) ::  harvest   !<
269
270      REAL(wp) ::  rsq      !<
271      REAL(wp) ::  v1       !<
272      REAL(wp) ::  v2       !<
273      REAL(wp), SAVE ::  g  !<
274
275      LOGICAL, SAVE ::  gaus_stored = .FALSE. !<
276!
277!--   We have an extra deviate handy, so return it, and unset the flag.
278      IF (gaus_stored)  THEN
279         harvest = g
280         gaus_stored = .FALSE.
281!
282!--   We don’t have an extra deviate handy, so pick two uniform numbers in the
283!--   square extending from -1 to +1 in each direction
284      ELSE
285         DO
286            CALL ran0_s(v1)
287            CALL ran0_s(v2)
288            v1 = 2.0_wp * v1 - 1.0_wp
289            v2 = 2.0_wp * v2 - 1.0_wp
290!
291!--         see if they are in the unit circle
292            rsq = v1**2 + v2**2
293!
294!--         otherwise try again.
295            IF (rsq > 0.0  .AND.  rsq < 1.0) EXIT
296         ENDDO
297!
298!--      Now make the Box-Muller transformation to get two normal deviates.
299!--      Return one and save the other for next time. Set flag.
300         rsq = SQRT(-2.0_sp * LOG(rsq)/rsq)
301         harvest = v1 * rsq
302         g = v2 * rsq
303         gaus_stored = .TRUE.
304      ENDIF
305
306   END SUBROUTINE gasdev_s
307
308!------------------------------------------------------------------------------!
309! Description:
310! ------------
311!> Initialize or reinitialize the random generator state space to vectors of size length.
312!> The saved variable seq is hashed (via calls to the module routine ran_hash)
313!> to create unique starting seeds, different for each vector component.
314!------------------------------------------------------------------------------!
315   SUBROUTINE ran_init( length )
316   
317      USE kinds
318     
319      IMPLICIT NONE
320     
321      INTEGER(isp), INTENT(IN) ::  length   !<
322   
323      INTEGER(isp) ::  hg    !<
324      INTEGER(isp) ::  hgm   !<
325      INTEGER(isp) ::  hgng  !<
326     
327      INTEGER(isp) ::  new   !<
328      INTEGER(isp) ::  j     !<
329      INTEGER(isp) ::  hgt   !<
330     
331      IF ( length < lenran ) RETURN !- Simply return if enough space is already allocated.
332     
333      hg=huge(1_isp)
334      hgm=-hg
335      hgng=hgm-1
336      hgt = hg
337     
338      !- The following lines check that kind value isp is in fact a 32-bit integer
339      !- with the usual properties that we expect it to have (under negation and wrap-around addition).
340      !- If all of these tests are satisfied, then the routines that use this module are portable,
341      !- even though they go beyond Fortran 90's integer model.
342     
343      IF  ( hg /= 2147483647 ) CALL ran_error('arithmetic assumption 1 failed')
344      IF  ( hgng >= 0 )        CALL ran_error('arithmetic assumption 2 failed')
345      IF  ( hgt+1 /= hgng )    CALL ran_error('arithmetic assumption 3 failed')
346      IF  ( not(hg) >= 0 )     CALL ran_error('arithmetic assumption 4 failed')
347      IF  ( not(hgng) < 0 )    CALL ran_error('arithmetic assumption 5 failed')
348      IF  ( hg+hgng >= 0 )     CALL ran_error('arithmetic assumption 6 failed')
349      IF  ( not(-1_isp) < 0 )  CALL ran_error('arithmetic assumption 7 failed')
350      IF  ( not(0_isp) >= 0 )  CALL ran_error('arithmetic assumption 8 failed')
351      IF  ( not(1_isp) >= 0 )  CALL ran_error('arithmetic assumption 9 failed')
352     
353      IF  ( lenran > 0) THEN                          !- Reallocate space, or ...
354     
355         ranseeds => reallocate( ranseeds, length, 5)
356         ranv => reallocate( ranv, length-1)
357         new = lenran+1
358         
359      ELSE                                            !- allocate space.
360     
361         ALLOCATE(ranseeds(length,5))
362         ALLOCATE(ranv(length-1))
363         new = 1   !- Index of first location not yet initialized.
364         amm = nearest(1.0_wp,-1.0_wp)/hgng
365         !- Use of nearest is to ensure that returned random deviates are strictly lessthan 1.0.
366         IF  (amm*hgng >= 1.0 .or. amm*hgng <= 0.0)                            &
367            CALL ran_error('arithmetic assumption 10 failed')
368           
369      END IF 
370     
371      !- Set starting values, unique by seq and vector component.
372      ranseeds(new:,1) = seq
373      ranseeds(new:,2:5)=spread(arth(new,1,size(ranseeds(new:,1))),2,4)
374     
375      DO j=1,4   !- Hash them.
376         CALL ran_hash(ranseeds(new:,j), ranseeds(new:,j+1))
377      END DO
378     
379      WHERE (ranseeds (new: ,1:3) < 0)                                         & 
380         ranseeds(new: ,1:3)=not(ranseeds(new: ,1:3))  !- Enforce nonnegativity.
381         
382      WHERE (ranseeds(new: ,4:5) == 0) ranseeds(new: ,4:5)=1 !- Enforce nonzero.
383     
384      IF  (new == 1) THEN !- Set scalar seeds.
385     
386         iran0 = ranseeds(1,1)
387         jran0 = ranseeds(1,2)
388         kran0 = ranseeds(1,3)
389         mran0 = ranseeds(1,4)
390         nran0 = ranseeds(1,5)
391         rans  = nran0
392         
393      END IF
394     
395      IF  (length > 1) THEN   !- Point to vector seeds.
396     
397         iran => ranseeds(2:,1)
398         jran => ranseeds(2:,2)
399         kran => ranseeds(2:,3)
400         mran => ranseeds(2:,4)
401         nran => ranseeds(2:,5)
402         ranv = nran
403         
404      END IF
405     
406      lenran = length
407     
408   END SUBROUTINE ran_init
409
410!------------------------------------------------------------------------------!
411! Description:
412! ------------
413!> User interface to release the workspace used by the random number routines.
414!------------------------------------------------------------------------------!
415   SUBROUTINE ran_deallocate
416   
417      IF  ( lenran > 0 ) THEN
418     
419         DEALLOCATE(ranseeds, ranv)
420         NULLIFY(ranseeds, ranv, iran, jran, kran, mran, nran)
421         lenran = 0
422         
423      END IF
424     
425   END SUBROUTINE ran_deallocate
426
427!------------------------------------------------------------------------------!
428! Description:
429! ------------
430!> User interface for seeding the random number routines.
431!> Syntax is exactly like Fortran 90's random_seed routine,
432!> with one additional argument keyword: random_sequence, set to any integer
433!> value, causes an immediate new initialization, seeded by that integer.
434!------------------------------------------------------------------------------!
435   SUBROUTINE random_seed_parallel( random_sequence, state_size, put, get )
436   
437      IMPLICIT NONE
438     
439      INTEGER(isp), OPTIONAL, INTENT(IN)  ::  random_sequence   !<
440      INTEGER(isp), OPTIONAL, INTENT(OUT) ::  state_size        !<
441     
442      INTEGER(isp), DIMENSION(:), OPTIONAL, INTENT(IN)  ::  put   !<
443      INTEGER(isp), DIMENSION(:), OPTIONAL, INTENT(OUT) ::  get   !<
444     
445      IF  ( present(state_size) ) THEN
446     
447         state_size = 5 * lenran
448         
449      ELSE IF  ( present(put) ) THEN
450     
451         IF  ( lenran == 0 ) RETURN
452         
453         ranseeds = reshape( put,shape(ranseeds) )
454         
455         WHERE (ranseeds(:,1:3) < 0) ranseeds(: ,1:3) = not(ranseeds(: ,1:3))
456         !- Enforce nonnegativity and nonzero conditions on any user-supplied seeds.
457         
458         WHERE (ranseeds(:,4:5) == 0) ranseeds(:,4:5) = 1
459         
460         iran0 = ranseeds(1,1)
461         jran0 = ranseeds(1,2)
462         kran0 = ranseeds(1,3)
463         mran0 = ranseeds(1,4)
464         nran0 = ranseeds(1,5)
465         
466      ELSE IF  ( present(get) ) THEN
467     
468         IF  (lenran == 0) RETURN
469         
470         ranseeds(1,1:5) = (/ iran0,jran0,kran0,mran0,nran0 /)
471         get = reshape( ranseeds, shape(get) )
472         
473      ELSE IF  ( present(random_sequence) ) THEN
474     
475         CALL ran_deallocate
476         seq = random_sequence
477         
478      END IF
479     
480   END SUBROUTINE random_seed_parallel
481
482!------------------------------------------------------------------------------!
483! Description:
484! ------------
485!> DES-like hashing of two 32-bit integers, using shifts,
486!> xor's, and adds to make the internal nonlinear function.
487!------------------------------------------------------------------------------!
488   SUBROUTINE ran_hash_v( il, ir )
489   
490      IMPLICIT NONE
491     
492      INTEGER(isp), DIMENSION(:), INTENT(INOUT) ::  il   !<
493      INTEGER(isp), DIMENSION(:), INTENT(INOUT) ::  ir   !<
494     
495      INTEGER(isp), DIMENSION(size(il)) ::  is   !<
496     
497      INTEGER(isp) :: j   !<
498     
499      DO j=1,4
500     
501         is = ir
502         ir = ieor( ir, ishft(ir,5) ) + 1422217823
503         ir = ieor( ir, ishft(ir,-16) ) + 1842055030
504         ir = ieor( ir, ishft(ir,9) ) + 80567781
505         ir = ieor( il, ir )
506         il = is
507      END DO
508     
509   END SUBROUTINE ran_hash_v
510
511!------------------------------------------------------------------------------!
512! Description:
513! ------------
514!> User interface to process error-messages
515!> produced by the random_number_parallel module
516!------------------------------------------------------------------------------!
517   SUBROUTINE ran_error(string)
518
519      USE control_parameters,                                                &
520        ONLY: message_string
521
522      CHARACTER(LEN=*), INTENT(IN) ::  string   !< Error message string
523
524      message_string = 'incompatible integer arithmetic: ' // TRIM( string )
525      CALL message( 'ran_init', 'PA0453', 1, 2, 0, 6, 0 )
526
527   END SUBROUTINE ran_error
528
529!------------------------------------------------------------------------------!
530! Description:
531! ------------
532!> Reallocates the generators state space "ranseeds" to vectors of size length.
533!------------------------------------------------------------------------------!
534   FUNCTION reallocate_iv( p, n )
535
536      USE control_parameters,                                                &
537        ONLY: message_string
538   
539      INTEGER(isp), DIMENSION(:), POINTER ::  p               !<
540      INTEGER(isp), DIMENSION(:), POINTER ::  reallocate_iv   !<
541     
542      INTEGER(isp), INTENT(IN) ::  n   !<
543     
544      INTEGER(isp) ::  nold   !<
545      INTEGER(isp) ::  ierr   !<
546     
547      ALLOCATE(reallocate_iv(n),stat=ierr)
548     
549      IF (ierr /= 0) THEN
550         message_string = 'problem in attempt to allocate memory'
551         CALL message( 'reallocate_iv', 'PA0454', 1, 2, 0, 6, 0 )
552      END IF
553
554      IF (.not. associated(p)) RETURN
555     
556      nold = size(p)
557     
558      reallocate_iv(1:min(nold,n)) = p(1:min(nold,n))
559     
560      DEALLOCATE(p)
561     
562   END FUNCTION reallocate_iv
563   
564   FUNCTION reallocate_im( p, n, m )
565
566      USE control_parameters,                                                &
567        ONLY: message_string
568   
569      INTEGER(isp), DIMENSION(:,:), POINTER ::  p               !<
570      INTEGER(isp), DIMENSION(:,:), POINTER ::  reallocate_im   !<
571     
572      INTEGER(isp), INTENT(IN) ::  m   !<
573      INTEGER(isp), INTENT(IN) ::  n   !<
574     
575      INTEGER(isp) ::  mold   !<
576      INTEGER(isp) ::  nold   !<
577      INTEGER(isp) ::  ierr   !<
578     
579      ALLOCATE(reallocate_im(n,m),stat=ierr)
580     
581      IF (ierr /= 0) THEN
582         message_string = 'problem in attempt to allocate memory'
583         CALL message( 'reallocate_im', 'PA0454', 1, 2, 0, 6, 0 )
584      END IF
585
586      IF (.not. associated(p)) RETURN
587     
588      nold = size(p,1)
589      mold = size(p,2)
590     
591      reallocate_im(1:min(nold,n),1:min(mold,m)) =                             &
592         p(1:min(nold,n),1:min(mold,m))
593         
594      DEALLOCATE(p)
595     
596   END FUNCTION reallocate_im
597
598!------------------------------------------------------------------------------!
599! Description:
600! ------------
601!> Reallocates the generators state space "ranseeds" to vectors of size length.
602!------------------------------------------------------------------------------!
603   FUNCTION arth_i(first,increment,n)
604   
605      INTEGER(isp), INTENT(IN) ::  first       !<
606      INTEGER(isp), INTENT(IN) ::  increment   !<
607      INTEGER(isp), INTENT(IN) ::  n           !<
608     
609      INTEGER(isp), DIMENSION(n) ::  arth_i    !<
610     
611      INTEGER(isp) ::  k      !<
612      INTEGER(isp) ::  k2     !<
613      INTEGER(isp) ::  temp   !<
614     
615      INTEGER(isp), PARAMETER ::  npar_arth=16   !<
616      INTEGER(isp), PARAMETER ::  npar2_arth=8   !<
617     
618      IF (n > 0) arth_i(1) = first
619     
620      IF (n <= npar_arth) THEN
621     
622         DO k=2,n
623            arth_i(k) = arth_i(k-1)+increment
624         END DO
625         
626      ELSE
627     
628         DO k=2,npar2_arth
629            arth_i(k) = arth_i(k-1) + increment
630         END DO
631         
632         temp = increment * npar2_arth
633         k = npar2_arth
634         
635         DO
636            IF (k >= n) EXIT
637            k2 = k + k
638            arth_i(k+1:min(k2,n)) = temp + arth_i(1:min(k,n-k))
639            temp = temp + temp
640            k = k2
641         END DO
642         
643      END IF
644     
645   END FUNCTION arth_i
646
647END MODULE random_generator_parallel
Note: See TracBrowser for help on using the repository browser.