Changeset 4592 for palm/trunk/SOURCE/random_generator_parallel_mod.f90
 Timestamp:
 Jul 8, 2020 9:56:29 AM (17 months ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

palm/trunk/SOURCE/random_generator_parallel_mod.f90
r4545 r4592 1 1 !> @file random_generator_parallel_mod.f90 2 ! !2 !! 3 3 ! This file is part of the PALM model system. 4 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/>. 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/>. 16 15 ! 17 16 ! Copyright 19972020 Leibniz Universitaet Hannover 18 !! 17 !! 18 ! 19 19 ! 20 20 ! Current revisions: 21 21 !  22 ! 23 ! 22 ! 23 ! 24 24 ! Former revisions: 25 25 !  26 26 ! $Id$ 27 ! Add generator (using parallel mode) returning gaussian distributed random 28 ! number 29 ! 27 ! File reformatted to follow the PALM coding standard 28 ! 29 ! 4545 20200522 13:17:57Z schwenkel 30 ! Add generator (using parallel mode) returning gaussian distributed random number 31 ! 30 32 ! 4438 20200303 20:49:28Z raasch 31 33 !  Rename variables to avoid confusion with subdomain grid indices 32 34 !  Some formatting adjustments 33 !  New routine to initialize spatial 1Darrays independent on the 2D random 34 ! number array 35 ! 35 !  New routine to initialize spatial 1Darrays independent on the 2D randomnumber array 36 ! 36 37 ! 4360 20200107 11:25:50Z suehring 37 38 ! Corrected "Former revisions" section 38 ! 39 ! 39 40 ! 3655 20190107 16:51:22Z knoop 40 ! unused variable removed41 ! Unused variable removed 41 42 ! 42 43 ! 1400 20140509 14:03:54Z knoop 43 44 ! Initial revision 44 ! 45 ! 45 ! 46 ! 47 !! 46 48 ! Description: 47 49 !  48 50 !> 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 32bit 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 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 32bit 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. 60 59 !> 61 60 !> This routine is taken from the "numerical recipies vol. 2" 62 !! 63 MODULE 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 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 128 131 129 132 CONTAINS 130 133 131 ! !134 !! 132 135 ! Description: 133 136 !  134 137 !> Initialize the parallel random number generator for a 1dimensional 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 ydirection 185 INTEGER(isp), INTENT(IN) :: nyn_l !< local upper subdomain bound index in ydirection 186 INTEGER(isp), INTENT(IN) :: nxl_l !< local lower subdomain bound index in xdirection 187 INTEGER(isp), INTENT(IN) :: nxr_l !< local upper subdomain bound index in xdirection 188 189 INTEGER(iwp) :: i !< grid index xdirection 190 INTEGER(iwp) :: j !< grid index ydirection 191 192 ! Allocate IDarray and statespacearray 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. 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 ydirection 185 INTEGER(isp), INTENT(IN) :: nyn_l !< local upper subdomain bound index in ydirection 186 INTEGER(isp), INTENT(IN) :: nxl_l !< local lower subdomain bound index in xdirection 187 INTEGER(isp), INTENT(IN) :: nxr_l !< local upper subdomain bound index in xdirection 188 189 INTEGER(iwp) :: i !< grid index xdirection 190 INTEGER(iwp) :: j !< grid index ydirection 191 192 ! Allocate IDarray and statespacearray 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. 226 224 !> The period of this generator is about 2.0 x 10^28, and it fully vectorizes. 227 225 !> 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 BoxMuller 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=hgm1 336 hgt = hg 337 338 ! The following lines check that kind value isp is in fact a 32bit integer 339 ! with the usual properties that we expect it to have (under negation and wraparound 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, length1) 357 new = lenran+1 358 359 ELSE ! allocate space. 360 361 ALLOCATE(ranseeds(length,5)) 362 ALLOCATE(ranv(length1)) 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 !! 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 BoxMuller 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 32bit integer with the usual 341 ! properties that we expect it to have (under negation and wraparound 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(length1) ) 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 !! 411 414 ! Description: 412 415 !  413 416 !> 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 usersupplied 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 !> DESlike hashing of two 32bit 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 errormessages 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(k1)+increment 624 END DO 625 626 ELSE 627 628 DO k=2,npar2_arth 629 arth_i(k) = arth_i(k1) + 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,nk)) 639 temp = temp + temp 640 k = k2 641 END DO 642 643 END IF 644 645 END FUNCTION arth_i 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 usersupplied 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 !> DESlike hashing of two 32bit 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 errormessages 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(k1) + increment 625 END DO 626 627 ELSE 628 629 DO k = 2, npar2_arth 630 arth_i(k) = arth_i(k1) + 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 646 647 647 648 END MODULE random_generator_parallel
Note: See TracChangeset
for help on using the changeset viewer.