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

Last change on this file since 4888 was 4828, checked in by Giersch, 3 years ago

Copyright updated to year 2021, interface pmc_sort removed to accelarate the nesting code

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