!!******************************************************************************
!!
!! This file is part of the AMUN source code, a program to perform
!! Newtonian or relativistic magnetohydrodynamical simulations on uniform or
!! adaptive mesh.
!!
!! Copyright (C) 2008-2021 Grzegorz Kowal <grzegorz@amuncode.org>
!!
!! This program 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.
!!
!! This program 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 this program. If not, see <http://www.gnu.org/licenses/>.
!!
!!******************************************************************************
!!
!! module: RANDOM
!!
!! This module provides subroutines to random number generators.
!!
!! References:
!!
!! [1] David Blackman and Sebastiano Vigna,
!! "Scrambled linear pseudorandom number generators", 2019
!! http://vigna.di.unimi.it/ftp/papers/ScrambledLinear.pdf
!! [2] http://prng.di.unimi.it/
!!
!!******************************************************************************
!
module random
implicit none
! random generator type
!
character(len = 32), save :: gentype = "same"
! variables to store seeds for random generator
!
integer , save :: kp = 0, np = 0
integer , save :: nseeds, lseed
integer(kind=8) , save :: state
integer(kind=8), dimension(:,:), allocatable, save :: seeds
! integer to real conversion parameters
!
real(kind=8), parameter :: i64tor8 = 2.0d+00**(-53)
! by default everything is private
!
private
! declare public subroutines
!
public :: initialize_random, finalize_random
public :: nseeds, set_seeds, get_seeds, gentype
public :: randuni, randsym, randnorz
contains
!
!===============================================================================
!
! subroutine INITIALIZE_RANDOM:
! ----------------------------
!
! Subroutine initializes random number generator.
!
! Arguments:
!
! gtype - the generator type ('same' or 'random');
! nthreads - the number of threads;
! nthread - the thread number;
! nproc - the process number;
! status - the call status;
!
!===============================================================================
!
subroutine initialize_random(gtype, nthreads, nthread, nproc, status)
implicit none
character(len=*), intent(in) :: gtype
integer , intent(in) :: nthreads, nthread, nproc
integer , intent(out) :: status
integer :: i
integer(kind=8), dimension(4) :: s
!-------------------------------------------------------------------------------
!
status = 0
! set local copies of the thread and process numbers,
! determine the number of seeds
!
kp = nthread
np = nproc
nseeds = nthreads
lseed = nseeds - 1
! allocate seeds
!
allocate(seeds(4, 0:lseed), stat = status)
if (status == 0) then
! prepare seeds depending on the type of generator
!
write(gentype,'(a)') trim(gtype)
select case(gentype)
case('random')
state = 1234567890
s(1) = splitmix64()
s(2) = splitmix64()
s(3) = splitmix64()
s(4) = splitmix64()
do i = 1, nseeds * np
call jump(s(:))
end do
seeds(:,0) = s(:)
do i = 1, lseed
call jump(s(:))
seeds(:,i) = s(:)
end do
case default
state = 1234567890
seeds(1,:) = splitmix64()
seeds(2,:) = splitmix64()
seeds(3,:) = splitmix64()
seeds(4,:) = splitmix64()
end select
end if
!-------------------------------------------------------------------------------
!
end subroutine initialize_random
!
!===============================================================================
!
! subroutine FINALIZE_RANDOM:
! --------------------------
!
! Subroutine releases memory allocated by random number generator variables.
!
! Arguments:
!
! status - the call status;
!
!===============================================================================
!
subroutine finalize_random(status)
implicit none
integer, intent(out) :: status
!-------------------------------------------------------------------------------
!
status = 0
if (allocated(seeds)) deallocate(seeds, stat=status)
!-------------------------------------------------------------------------------
!
end subroutine finalize_random
!
!===============================================================================
!
! subroutine SET_SEEDS:
! --------------------
!
! Subroutine sets the seeds.
!
! Arguments:
!
! n - the number of seeds;
! sds - the new seeds;
! gen - the flag indication if the seeds have to be regenerated;
!
!===============================================================================
!
subroutine set_seeds(n, sds, gen)
implicit none
integer , intent(in) :: n
integer(kind=8), dimension(4,n), intent(in) :: sds
logical , intent(in) :: gen
integer :: i
integer(kind=8), dimension(4) :: s
!-------------------------------------------------------------------------------
!
! set the seeds only if the input array and seeds have the same sizes
!
if (gen .or. n /= nseeds) then
select case(gentype)
case('random')
state = 1234567890
s(1) = splitmix64()
s(2) = splitmix64()
s(3) = splitmix64()
s(4) = splitmix64()
do i = 1, nseeds * np
call jump(s(:))
end do
seeds(:,0) = s(:)
do i = 1, lseed
call jump(s(:))
seeds(:,i) = s(:)
end do
case default
state = 1234567890
seeds(1,:) = splitmix64()
seeds(2,:) = splitmix64()
seeds(3,:) = splitmix64()
seeds(4,:) = splitmix64()
end select
else
seeds(:,:) = sds(:,:)
end if
!-------------------------------------------------------------------------------
!
end subroutine set_seeds
!
!===============================================================================
!
! subroutine GET_SEEDS:
! --------------------
!
! Subroutine returns the seeds.
!
! Arguments:
!
! sds - the return seeds array;
!
!===============================================================================
!
subroutine get_seeds(sds)
implicit none
integer(kind=8), dimension(4,nseeds), intent(out) :: sds
!-------------------------------------------------------------------------------
!
if (size(sds,1) == 4 .and. size(sds,2) == nseeds) then
sds(:,:) = seeds(:,:)
end if
!-------------------------------------------------------------------------------
!
end subroutine get_seeds
!
!===============================================================================
!
! function SPLITMIX64:
! -------------------
!
! Pseudorandom Number Generator for floating point numbers based on
! SplitMix64 algorithm [1], used for initialization of xoshiro256p()
! generator.
!
! References:
!
! [1] David Blackman and Sebastiano Vigna,
! "Scrambled linear pseudorandom number generators", 2019
! http://vigna.di.unimi.it/ftp/papers/ScrambledLinear.pdf
!
!===============================================================================
!
integer(kind=8) function splitmix64()
implicit none
integer(kind=8), parameter :: a = -7046029254386353131_8
integer(kind=8), parameter :: b = -4658895280553007687_8
integer(kind=8), parameter :: c = -7723592293110705685_8
!-------------------------------------------------------------------------------
!
state = state + a
splitmix64 = state
splitmix64 = ieor(splitmix64, ishft(splitmix64, -30)) * b
splitmix64 = ieor(splitmix64, ishft(splitmix64, -27)) * c
splitmix64 = ieor(splitmix64, ishft(splitmix64, -31))
return
!-------------------------------------------------------------------------------
!
end function splitmix64
!
!===============================================================================
!
! function XOSHIRO256P:
! --------------------
!
! Pseudorandom Number Generator for floating point numbers based on
! xoshiro256+ algorithm [1]. It return a 256-bit random integer.
!
! References:
!
! [1] David Blackman and Sebastiano Vigna,
! "Scrambled linear pseudorandom number generators", 2019
! http://vigna.di.unimi.it/ftp/papers/ScrambledLinear.pdf
!
!===============================================================================
!
integer(kind=8) function xoshiro256p()
implicit none
integer(kind=8) :: temp
!-------------------------------------------------------------------------------
!
xoshiro256p = ishft(seeds(1,kp) + seeds(4,kp), -11)
temp = ishft(seeds(2,kp), 17)
seeds(3,kp) = ieor(seeds(3,kp), seeds(1,kp))
seeds(4,kp) = ieor(seeds(4,kp), seeds(2,kp))
seeds(2,kp) = ieor(seeds(2,kp), seeds(3,kp))
seeds(1,kp) = ieor(seeds(1,kp), seeds(4,kp))
seeds(3,kp) = ieor(seeds(3,kp), temp)
seeds(4,kp) = ior(ishft(seeds(4,kp), 45), ishft(seeds(4,kp), -19))
return
!-------------------------------------------------------------------------------
!
end function xoshiro256p
!
!===============================================================================
!
! subroutine JUMP:
! ---------------
!
! The jump function for the xoshiro256p() generator used to generate
! non-overlapping seed subsequences for parallel computations.
!
! Arguments:
!
! sds - temporary seeds array;
!
! References:
!
! [1] David Blackman and Sebastiano Vigna,
! "Scrambled linear pseudorandom number generators", 2019
! http://vigna.di.unimi.it/ftp/papers/ScrambledLinear.pdf
!
!===============================================================================
!
subroutine jump(sds)
implicit none
integer(kind=8), dimension(4), intent(inout) :: sds
integer(kind=8), dimension(4) :: s
integer(kind=8) :: temp
integer :: i, b
integer(kind=8), dimension(4), parameter :: jmp = [ 1733541517147835066_8, &
-3051731464161248980_8, &
-6244198995065845334_8, &
4155657270789760540_8 ]
!-------------------------------------------------------------------------------
!
s(:) = 0_8
do i = 1, 4
do b = 1, 64
if (iand(jmp(i), ishft(1_8, b)) /= 0_8) then
s(1) = ieor(s(1), sds(1))
s(2) = ieor(s(2), sds(2))
s(3) = ieor(s(3), sds(3))
s(4) = ieor(s(4), sds(4))
end if
temp = ishft(sds(2), 17)
sds(3) = ieor(sds(3), sds(1))
sds(4) = ieor(sds(4), sds(2))
sds(2) = ieor(sds(2), sds(3))
sds(1) = ieor(sds(1), sds(4))
sds(3) = ieor(sds(3), temp)
sds(4) = ior(ishft(sds(4), 45), ishft(sds(4), -19))
end do
end do
sds(:) = s(:)
return
!-------------------------------------------------------------------------------
!
end subroutine jump
!
!===============================================================================
!
! function RANDUNI:
! ----------------
!
! Function generates uniformly distributed floating point real random number
! in range 0.0..1.0 using xoshiro256p() generator.
!
!===============================================================================
!
real(kind=8) function randuni()
implicit none
!-------------------------------------------------------------------------------
!
randuni = i64tor8 * xoshiro256p()
return
!-------------------------------------------------------------------------------
!
end function randuni
!
!===============================================================================
!
! function RANDSYM:
! ----------------
!
! Function generates uniformly distributed floating point random number
! in range -1.0..1.0 using xoshiro256p() generator.
!
!===============================================================================
!
real(kind=8) function randsym()
implicit none
!-------------------------------------------------------------------------------
!
randsym = 2.0d+00 * (i64tor8 * xoshiro256p() - 5.0d-01)
return
!-------------------------------------------------------------------------------
!
end function randsym
!
!===============================================================================
!
! function RANDNORZ:
! -----------------
!
! Function generates complex floating point number with a normal
! distribution (μ = 0, σ = 1) using Marsaglia polar method [1].
!
! References:
!
! [1] Marsaglia, G. & Bray, T. A.,
! "A Convenient Method for Generating Normal Variables",
! SIAM Review, vol. 6, no. 3, pp. 260-264,
! https://doi.org/10.1137/1006063
!
!===============================================================================
!
complex(kind=8) function randnorz()
implicit none
logical :: c
real(kind=8) :: u, v, s
!-------------------------------------------------------------------------------
!
c = .true.
do while(c)
u = randsym()
v = randsym()
s = u**2 + v**2
c = s <= 0.0d+00 .or. s >= 1.0d+00
end do
randnorz = sqrt(- 2.0d+00 * log(s) / s) * cmplx(u, v, kind=8)
return
!-------------------------------------------------------------------------------
!
end function randnorz
!===============================================================================
!
end module random