ran.fpp Source File


Contents

Source Code


Source Code

# 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).
  function ranf (seed)
# if RANDOM == _RANMT_
    use mt19937, only: sgrnd, grnd
# endif
    implicit none
    integer, intent(in), optional :: seed
    real :: ranf
    integer :: l
    integer, allocatable :: seed_in(:)
# 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
    integer :: i, nseed, dt(8)
    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

  !> 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
# if RANDOM == _RANMT_
    allocate(seed(1))
    seed = initial_seed
# else
    ! 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
# endif
  end function create_seed_from_single_integer

  !> 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