# include "define.inc" !> A wrapper module for random number generator. !> This module provides real function ranf !> using intrinsic random_number/random_seed or !> Mersenne Twister 19937 (see [[mt19937.f90]]). module ran implicit none private public :: ranf, set_seed_from_single_integer public :: get_rnd_seed_length,get_rnd_seed,init_ranf contains !> Returns a pseudorandom number range in [0., 1.), (i.e. 0<= ranf < 1). !> The generator is initialized with the given seed if passed, !> otherwise uses the seed already set (default or otherwise). real function ranf (seed) # if RANDOM == _RANMT_ use mt19937, only: sgrnd, grnd # endif implicit none integer, intent(in), optional :: seed # if RANDOM != _RANMT_ integer :: l integer, allocatable :: seed_in(:) # endif # if RANDOM == _RANMT_ if (present(seed)) call sgrnd(seed) ranf = grnd() # else if (present(seed)) then call random_seed(size=l) allocate(seed_in(l)) ! @note This is probably not ideal, instead we could use ! [[set_seed_from_single_integer]] instead seed_in(:)=seed call random_seed(put=seed_in) endif call random_number(ranf) # endif end function ranf !> Gets the length of the integer vector for the random number !> generator seed. This is always 1 for the Mersenne case but !> otherwise depends on the compiler version. function get_rnd_seed_length () result (l) implicit none integer :: l # if RANDOM == _RANMT_ l=1 # else call random_seed(size=l) # endif end function get_rnd_seed_length !> Returns the current value of the generator seed. !> This is not currently supported by the Mersenne !> generator so in this case we return 0. subroutine get_rnd_seed(seed) implicit none integer, dimension(:), intent(out) :: seed # if RANDOM == _RANMT_ !NOTE: If MT is generator, more coding needs to be done to ! extract seed, so this should be run using only randomize seed=0 # else call random_seed(get=seed) # endif end subroutine get_rnd_seed !> Seed the choosen random number generator. !> If randomize=T, a random seed based on date and time is used. !> (this randomizing procedure is proposed in GNU gfortran manual.) !> Otherwise, it sets the seed using init_seed !> In either case, it sets `init_seed` to the initial seed used on return. subroutine init_ranf(randomize, init_seed) use constants, only: kind_id # if RANDOM == _RANMT_ use mt19937, only: sgrnd # endif implicit none logical, intent(in) :: randomize integer, intent(in out), dimension(:) :: init_seed # if RANDOM != _RANMT_ integer :: i, nseed # endif integer, dimension(8) :: dt integer (kind=kind_id) :: t if (randomize) then call system_clock(t) if (t == 0) then call date_and_time(values=dt) t = (dt(1) - 1970) * 365_kind_id * 24 * 60 * 60 * 1000 & + dt(2) * 31_kind_id * 24 * 60 * 60 * 1000 & + dt(3) * 24_kind_id * 60 * 60 * 1000 & + dt(5) * 60 * 60 * 1000 & + dt(6) * 60 * 1000 + dt(7) * 1000 & + dt(8) end if else t = 0 end if # if RANDOM == _RANMT_ if (randomize) then init_seed(1) = lcg(t) endif call sgrnd(init_seed(1)) # else if (randomize) then ! We should probably check that size(init_seed) >= nseed call random_seed(size=nseed) do i = 1, nseed init_seed(i) = lcg(t) end do endif call random_seed(put=init_seed) # endif end subroutine init_ranf # if RANDOM != _RANMT_ !> Returns an array of integers to use as a seed to the !> random number generator given a single integer as input. !> !> Currently we just use the [[lcg]] generator to generate !> a series of values to use for the seed. We could also !> set the actual random seed from this and then use this !> to generate a final series of seeds. This is probably not !> really necessary. !> !> For our mt19937 implementation we just set the seed !> equal to initial_seed function create_seed_from_single_integer(initial_seed) result(seed) use constants, only: kind_id implicit none integer, intent(in) :: initial_seed integer, dimension(:), allocatable :: seed integer :: seed_size, i integer(kind = kind_id) :: current ! Get the required state size for the seed call random_seed(size = seed_size) allocate(seed(seed_size)) ! Use lcg to generate a series of random numbers given ! the initial seed current = initial_seed do i = 1, seed_size seed(i) = lcg(current) current = seed(i) end do end function create_seed_from_single_integer # endif !> Sets the seed for the random number generator provided by the !> standard based on a single integer. We use !> [[create_seed_from_single_integer]] to ensure we have enough values !> in the series rather than just duplicating this number. subroutine set_seed_from_single_integer(initial_seed) # if RANDOM == _RANMT_ use mt19937, only: sgrnd # endif implicit none integer, intent(in) :: initial_seed # if RANDOM == _RANMT_ call sgrnd(initial_seed) # else call random_seed(put = create_seed_from_single_integer(initial_seed)) # endif end subroutine set_seed_from_single_integer !> Simple PRNG derived from [gfortran !> documentation](https://gcc.gnu.org/onlinedocs/gcc-6.4.0/gfortran/RANDOM_005fSEED.html) !> Intended to be used to seed a better PRNG. Probably shouldn't be used !> directly (hence why this is not public). !> !> @warning (?) the sign of the output will be the same as the input !> and hence has a restricted range, perhaps we should consider !> subtracting sign(s)*huge(s)/2 ? function lcg(s) use constants, only: kind_id, kind_is implicit none integer :: lcg integer (kind=kind_id), intent(in out) :: s if (kind_id > 0) then if (s == 0) then s = 104729 else s = mod(s, 4294967296_kind_id) end if s = mod(s * 279470273_kind_id, 4294967291_kind_id) lcg = int(mod(s, int(huge(0), kind_id)), kind(0)) else lcg = mod(16807*int(s,kind_is), 2147483647) end if end function lcg end module ran