!> @file random_generator_parallel_mod.f90 !------------------------------------------------------------------------------! ! This file is part of PALM. ! ! PALM is free software: you can redistribute it and/or modify it under the ! terms of the GNU General Public License as published by the Free Software ! Foundation, either version 3 of the License, or (at your option) any later ! version. ! ! PALM is distributed in the hope that it will be useful, but WITHOUT ANY ! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR ! A PARTICULAR PURPOSE. See the GNU General Public License for more details. ! ! You should have received a copy of the GNU General Public License along with ! PALM. If not, see . ! ! Copyright 1997-2017 Leibniz Universitaet Hannover !------------------------------------------------------------------------------! ! ! Current revisions: ! ----------------- ! ! ! Former revisions: ! ----------------- ! $Id: random_generator_parallel_mod.f90 2255 2017-06-07 11:58:09Z sward $ ! Implemented PALM specific error message handling ! ! 2173 2017-03-08 15:56:54Z knoop ! ! 2172 2017-03-08 15:55:25Z knoop ! Bugfix: added global initialization routine and removed global array ! ! 2144 2017-02-07 15:32:45Z knoop ! Bugfix: compile time parameter replacement prevented ! ! 2000 2016-08-20 18:09:15Z knoop ! Forced header and separation lines into 80 columns ! ! 1850 2016-04-08 13:29:27Z maronga ! Module renamed ! ! ! 1682 2015-10-07 23:56:08Z knoop ! Code annotations made doxygen readable ! ! 1400 2014-05-09 14:03:54Z knoop ! Initial revision ! ! ! Description: ! ------------ !> This module contains and supports the random number generating routine ran_parallel. !> ran_parallel returns a uniform random deviate between 0.0 and 1.0 !> (exclusive of the end point values). !> Additionally it provides the generator with five integer for use as initial state space. !> The first tree integers (iran, jran, kran) are maintained as non negative values, !> while the last two (mran, nran) have 32-bit nonzero values. !> Also provided by this module is support for initializing or reinitializing !> the state space to a desired standard sequence number, hashing the initial !> values to random values, and allocating and deallocating the internal workspace !> Random number generator, produces numbers equally distributed in interval !> !> This routine is taken from the "numerical recipies vol. 2" !------------------------------------------------------------------------------! MODULE random_generator_parallel USE kinds IMPLICIT NONE PRIVATE PUBLIC random_number_parallel, random_seed_parallel, random_dummy, & id_random_array, seq_random_array, init_parallel_random_generator INTEGER(isp), SAVE :: lenran=0 !< INTEGER(isp), SAVE :: seq=0 !< INTEGER(isp), SAVE :: iran0 !< INTEGER(isp), SAVE :: jran0 !< INTEGER(isp), SAVE :: kran0 !< INTEGER(isp), SAVE :: mran0 !< INTEGER(isp), SAVE :: nran0 !< INTEGER(isp), SAVE :: rans !< INTEGER(isp), DIMENSION(:, :), POINTER, SAVE :: ranseeds !< INTEGER(isp), DIMENSION(:), POINTER, SAVE :: iran !< INTEGER(isp), DIMENSION(:), POINTER, SAVE :: jran !< INTEGER(isp), DIMENSION(:), POINTER, SAVE :: kran !< INTEGER(isp), DIMENSION(:), POINTER, SAVE :: mran !< INTEGER(isp), DIMENSION(:), POINTER, SAVE :: nran !< INTEGER(isp), DIMENSION(:), POINTER, SAVE :: ranv !< INTEGER(isp), DIMENSION(:,:), ALLOCATABLE :: id_random_array !< INTEGER(isp), DIMENSION(:,:,:), ALLOCATABLE :: seq_random_array !< REAL(wp), SAVE :: amm !< REAL(wp) :: random_dummy=0.0 !< INTERFACE init_parallel_random_generator MODULE PROCEDURE init_parallel_random_generator END INTERFACE INTERFACE random_number_parallel MODULE PROCEDURE ran0_s END INTERFACE INTERFACE random_seed_parallel MODULE PROCEDURE random_seed_parallel END INTERFACE INTERFACE ran_hash MODULE PROCEDURE ran_hash_v END INTERFACE INTERFACE reallocate MODULE PROCEDURE reallocate_iv,reallocate_im END INTERFACE INTERFACE arth MODULE PROCEDURE arth_i END INTERFACE CONTAINS !------------------------------------------------------------------------------! ! Description: ! ------------ !> Initialize the parallel random number generator for a specific subdomain !------------------------------------------------------------------------------! SUBROUTINE init_parallel_random_generator(nx, ny, nys, nyn, nxl, nxr) USE kinds USE control_parameters, & ONLY: ensemble_member_nr IMPLICIT NONE INTEGER(isp), INTENT(IN) :: nx !< INTEGER(isp), INTENT(IN) :: ny !< INTEGER(isp), INTENT(IN) :: nys !< INTEGER(isp), INTENT(IN) :: nyn !< INTEGER(isp), INTENT(IN) :: nxl !< INTEGER(isp), INTENT(IN) :: nxr !< INTEGER(iwp) :: i !< INTEGER(iwp) :: j !< !-- Allocate ID-array and state-space-array ALLOCATE ( seq_random_array(5,nys:nyn,nxl:nxr) ) ALLOCATE ( id_random_array(nys:nyn,nxl:nxr) ) seq_random_array = 0 id_random_array = 0 !-- Asigning an ID to every vertical gridpoint column !-- dependig on the ensemble run number. DO i=nxl, nxr DO j=nys, nyn id_random_array(j,i) = j*(nx+1.0_wp) + i + 1.0_wp + 1E6 * & ( ensemble_member_nr - 1000 ) ENDDO ENDDO !-- Initializing with random_seed_parallel for every vertical !-- gridpoint column. random_dummy=0 DO i = nxl, nxr DO j = nys, nyn CALL random_seed_parallel (random_sequence=id_random_array(j, i)) CALL random_number_parallel (random_dummy) CALL random_seed_parallel (get=seq_random_array(:, j, i)) ENDDO ENDDO END SUBROUTINE init_parallel_random_generator !------------------------------------------------------------------------------! ! Description: ! ------------ !> Lagged Fibonacci generator combined with a Marsaglia shiftsequence. !> Returns as harvest a uniform random deviate between 0.0 and 1.0 (exclusive of the end point values). !> This generator has the same calling and initialization conventions as Fortran 90's random_number routine. !> Use random_seed_parallel to initialize or reinitialize to a particular sequence. !> The period of this generator is about 2.0 x 10^28, and it fully vectorizes. !> Validity of the integer model assumed by this generator is tested at initialization. !------------------------------------------------------------------------------! SUBROUTINE ran0_s(harvest) USE kinds IMPLICIT NONE REAL(wp), INTENT(OUT) :: harvest !< IF (lenran < 1) CALL ran_init(1) !- Initialization routine in ran_state. !- Update Fibonacci generator, which has period p^2 + p + 1, p = 2^31 - 69. rans = iran0 - kran0 IF (rans < 0) rans = rans + 2147483579_isp iran0 = jran0 jran0 = kran0 kran0 = rans nran0 = ieor( nran0, ishft (nran0, 13) ) !- Update Marsaglia shift sequence with period 2^32 - 1. nran0 = ieor( nran0, ishft (nran0, -17) ) nran0 = ieor( nran0, ishft (nran0, 5) ) rans = ieor( nran0, rans ) !- Combine the generators. harvest = amm * merge( rans, not(rans), rans < 0 ) !- Make the result positive definite (note that amm is negative). END SUBROUTINE ran0_s !------------------------------------------------------------------------------! ! Description: ! ------------ !> Initialize or reinitialize the random generator state space to vectors of size length. !> The saved variable seq is hashed (via calls to the module routine ran_hash) !> to create unique starting seeds, different for each vector component. !------------------------------------------------------------------------------! SUBROUTINE ran_init( length ) USE kinds IMPLICIT NONE INTEGER(isp), INTENT(IN) :: length !< INTEGER(isp) :: hg !< INTEGER(isp) :: hgm !< INTEGER(isp) :: hgng !< INTEGER(isp) :: new !< INTEGER(isp) :: j !< INTEGER(isp) :: hgt !< IF ( length < lenran ) RETURN !- Simply return if enough space is already allocated. hg=huge(1_isp) hgm=-hg hgng=hgm-1 hgt = hg !- The following lines check that kind value isp is in fact a 32-bit integer !- with the usual properties that we expect it to have (under negation and wrap-around addition). !- If all of these tests are satisfied, then the routines that use this module are portable, !- even though they go beyond Fortran 90's integer model. IF ( hg /= 2147483647 ) CALL ran_error('arithmetic assumption 1 failed') IF ( hgng >= 0 ) CALL ran_error('arithmetic assumption 2 failed') IF ( hgt+1 /= hgng ) CALL ran_error('arithmetic assumption 3 failed') IF ( not(hg) >= 0 ) CALL ran_error('arithmetic assumption 4 failed') IF ( not(hgng) < 0 ) CALL ran_error('arithmetic assumption 5 failed') IF ( hg+hgng >= 0 ) CALL ran_error('arithmetic assumption 6 failed') IF ( not(-1_isp) < 0 ) CALL ran_error('arithmetic assumption 7 failed') IF ( not(0_isp) >= 0 ) CALL ran_error('arithmetic assumption 8 failed') IF ( not(1_isp) >= 0 ) CALL ran_error('arithmetic assumption 9 failed') IF ( lenran > 0) THEN !- Reallocate space, or ... ranseeds => reallocate( ranseeds, length, 5) ranv => reallocate( ranv, length-1) new = lenran+1 ELSE !- allocate space. ALLOCATE(ranseeds(length,5)) ALLOCATE(ranv(length-1)) new = 1 !- Index of first location not yet initialized. amm = nearest(1.0_wp,-1.0_wp)/hgng !- Use of nearest is to ensure that returned random deviates are strictly lessthan 1.0. IF (amm*hgng >= 1.0 .or. amm*hgng <= 0.0) & CALL ran_error('arithmetic assumption 10 failed') END IF !- Set starting values, unique by seq and vector component. ranseeds(new:,1) = seq ranseeds(new:,2:5)=spread(arth(new,1,size(ranseeds(new:,1))),2,4) DO j=1,4 !- Hash them. CALL ran_hash(ranseeds(new:,j), ranseeds(new:,j+1)) END DO WHERE (ranseeds (new: ,1:3) < 0) & ranseeds(new: ,1:3)=not(ranseeds(new: ,1:3)) !- Enforce nonnegativity. WHERE (ranseeds(new: ,4:5) == 0) ranseeds(new: ,4:5)=1 !- Enforce nonzero. IF (new == 1) THEN !- Set scalar seeds. iran0 = ranseeds(1,1) jran0 = ranseeds(1,2) kran0 = ranseeds(1,3) mran0 = ranseeds(1,4) nran0 = ranseeds(1,5) rans = nran0 END IF IF (length > 1) THEN !- Point to vector seeds. iran => ranseeds(2:,1) jran => ranseeds(2:,2) kran => ranseeds(2:,3) mran => ranseeds(2:,4) nran => ranseeds(2:,5) ranv = nran END IF lenran = length END SUBROUTINE ran_init !------------------------------------------------------------------------------! ! Description: ! ------------ !> User interface to release the workspace used by the random number routines. !------------------------------------------------------------------------------! SUBROUTINE ran_deallocate IF ( lenran > 0 ) THEN DEALLOCATE(ranseeds, ranv) NULLIFY(ranseeds, ranv, iran, jran, kran, mran, nran) lenran = 0 END IF END SUBROUTINE ran_deallocate !------------------------------------------------------------------------------! ! Description: ! ------------ !> User interface for seeding the random number routines. !> Syntax is exactly like Fortran 90's random_seed routine, !> with one additional argument keyword: random_sequence, set to any integer !> value, causes an immediate new initialization, seeded by that integer. !------------------------------------------------------------------------------! SUBROUTINE random_seed_parallel( random_sequence, state_size, put, get ) IMPLICIT NONE INTEGER(isp), OPTIONAL, INTENT(IN) :: random_sequence !< INTEGER(isp), OPTIONAL, INTENT(OUT) :: state_size !< INTEGER(isp), DIMENSION(:), OPTIONAL, INTENT(IN) :: put !< INTEGER(isp), DIMENSION(:), OPTIONAL, INTENT(OUT) :: get !< IF ( present(state_size) ) THEN state_size = 5 * lenran ELSE IF ( present(put) ) THEN IF ( lenran == 0 ) RETURN ranseeds = reshape( put,shape(ranseeds) ) WHERE (ranseeds(:,1:3) < 0) ranseeds(: ,1:3) = not(ranseeds(: ,1:3)) !- Enforce nonnegativity and nonzero conditions on any user-supplied seeds. WHERE (ranseeds(:,4:5) == 0) ranseeds(:,4:5) = 1 iran0 = ranseeds(1,1) jran0 = ranseeds(1,2) kran0 = ranseeds(1,3) mran0 = ranseeds(1,4) nran0 = ranseeds(1,5) ELSE IF ( present(get) ) THEN IF (lenran == 0) RETURN ranseeds(1,1:5) = (/ iran0,jran0,kran0,mran0,nran0 /) get = reshape( ranseeds, shape(get) ) ELSE IF ( present(random_sequence) ) THEN CALL ran_deallocate seq = random_sequence END IF END SUBROUTINE random_seed_parallel !------------------------------------------------------------------------------! ! Description: ! ------------ !> DES-like hashing of two 32-bit integers, using shifts, !> xor's, and adds to make the internal nonlinear function. !------------------------------------------------------------------------------! SUBROUTINE ran_hash_v( il, ir ) IMPLICIT NONE INTEGER(isp), DIMENSION(:), INTENT(INOUT) :: il !< INTEGER(isp), DIMENSION(:), INTENT(INOUT) :: ir !< INTEGER(isp), DIMENSION(size(il)) :: is !< INTEGER(isp) :: j !< DO j=1,4 is = ir ir = ieor( ir, ishft(ir,5) ) + 1422217823 ir = ieor( ir, ishft(ir,-16) ) + 1842055030 ir = ieor( ir, ishft(ir,9) ) + 80567781 ir = ieor( il, ir ) il = is END DO END SUBROUTINE ran_hash_v !------------------------------------------------------------------------------! ! Description: ! ------------ !> User interface to process error-messages !> produced by the random_number_parallel module !------------------------------------------------------------------------------! SUBROUTINE ran_error(string) USE control_parameters, & ONLY: message_string CHARACTER(LEN=*), INTENT(IN) :: string !< Error message string message_string = 'incompatible integer arithmetic: ' // TRIM( string ) CALL message( 'ran_init', 'PA0453', 1, 2, 0, 6, 0 ) END SUBROUTINE ran_error !------------------------------------------------------------------------------! ! Description: ! ------------ !> Reallocates the generators state space "ranseeds" to vectors of size length. !------------------------------------------------------------------------------! FUNCTION reallocate_iv( p, n ) USE control_parameters, & ONLY: message_string INTEGER(isp), DIMENSION(:), POINTER :: p !< INTEGER(isp), DIMENSION(:), POINTER :: reallocate_iv !< INTEGER(isp), INTENT(IN) :: n !< INTEGER(isp) :: nold !< INTEGER(isp) :: ierr !< ALLOCATE(reallocate_iv(n),stat=ierr) IF (ierr /= 0) THEN message_string = 'problem in attempt to allocate memory' CALL message( 'reallocate_iv', 'PA0454', 1, 2, 0, 6, 0 ) END IF IF (.not. associated(p)) RETURN nold = size(p) reallocate_iv(1:min(nold,n)) = p(1:min(nold,n)) DEALLOCATE(p) END FUNCTION reallocate_iv FUNCTION reallocate_im( p, n, m ) USE control_parameters, & ONLY: message_string INTEGER(isp), DIMENSION(:,:), POINTER :: p !< INTEGER(isp), DIMENSION(:,:), POINTER :: reallocate_im !< INTEGER(isp), INTENT(IN) :: m !< INTEGER(isp), INTENT(IN) :: n !< INTEGER(isp) :: mold !< INTEGER(isp) :: nold !< INTEGER(isp) :: ierr !< ALLOCATE(reallocate_im(n,m),stat=ierr) IF (ierr /= 0) THEN message_string = 'problem in attempt to allocate memory' CALL message( 'reallocate_im', 'PA0454', 1, 2, 0, 6, 0 ) END IF IF (.not. associated(p)) RETURN nold = size(p,1) mold = size(p,2) reallocate_im(1:min(nold,n),1:min(mold,m)) = & p(1:min(nold,n),1:min(mold,m)) DEALLOCATE(p) END FUNCTION reallocate_im !------------------------------------------------------------------------------! ! Description: ! ------------ !> Reallocates the generators state space "ranseeds" to vectors of size length. !------------------------------------------------------------------------------! FUNCTION arth_i(first,increment,n) INTEGER(isp), INTENT(IN) :: first !< INTEGER(isp), INTENT(IN) :: increment !< INTEGER(isp), INTENT(IN) :: n !< INTEGER(isp), DIMENSION(n) :: arth_i !< INTEGER(isp) :: k !< INTEGER(isp) :: k2 !< INTEGER(isp) :: temp !< INTEGER(isp), PARAMETER :: npar_arth=16 !< INTEGER(isp), PARAMETER :: npar2_arth=8 !< IF (n > 0) arth_i(1) = first IF (n <= npar_arth) THEN DO k=2,n arth_i(k) = arth_i(k-1)+increment END DO ELSE DO k=2,npar2_arth arth_i(k) = arth_i(k-1) + increment END DO temp = increment * npar2_arth k = npar2_arth DO IF (k >= n) EXIT k2 = k + k arth_i(k+1:min(k2,n)) = temp + arth_i(1:min(k,n-k)) temp = temp + temp k = k2 END DO END IF END FUNCTION arth_i END MODULE random_generator_parallel