Initializing RANDOM_NUMBER

3

To which the user says, “It’s your job to give me random numbers.” :slight_smile: But NAG has thought about this:

Question How should RANDOM_SEED be used to initialise the random number generator?

…

In the case of wanting the same sequence with every execution, you can use RANDOM_SEED with your own seed data. The Mersenne Twister state is a random bit sequence and it should therefore have approximately equal numbers of zero and one bits. If you do not have 630 32-bit integers worth of random data, you should set the remaining elements of the seed array to zero so that the generator will know that these do not contain any useful entropy.

My modified code that handles this case is

module random_mod
implicit none
contains
subroutine set_random_seed(fixed,offset,nburn_random)
logical, intent(in)           :: fixed        ! if .true., use a fixed seed
integer, intent(in), optional :: offset       ! constant added to generated seeds
integer, intent(in), optional :: nburn_random ! # of variates of RANDOM_NUMBER to burn off
integer, parameter            :: nburn = 100, nseeds_max = 8, seeds_vec(nseeds_max) = [675495868, 747035054, &
                                 914611756, 719350065, 933270595, 962123778, 509555460, 423028787] ! random integers
integer                       :: iburn,nseeds,itime,offset_,nburn_random_
integer, allocatable          :: seed(:)
real(kind=kind(1.0d0))        :: xran,xseed(nburn + nseeds_max)
nburn_random_ = 100
offset_ = 0
if (present(nburn_random)) nburn_random_ = nburn_random
if (present(offset)) offset_ = offset
call random_seed(size=nseeds)
allocate(seed(nseeds))
if (fixed) then
   if (nseeds > nseeds_max) then
      seed(:nseeds_max) = seeds_vec + offset_
      seed(nseeds_max + 1:) = 0 ! NAG says to set the remaining elements of the seed array to zero to indicate zero entropy
   end if
   seed = seeds_vec(:nseeds) + offset_
else ! set seed based on current time
   call system_clock(itime)
   seed = itime
   call random_seed(put=seed)
   call random_number(xseed)
   seed = xseed(nburn+1:) * 1000000000 + offset_
end if
call random_seed(put=seed)
do iburn=1,nburn_random_
   call random_number(xran)
end do
end subroutine set_random_seed
end module random_mod

program xset_random_seed
use random_mod, only: set_random_seed
implicit none
integer, parameter :: dp = kind(1.0d0)
integer, parameter :: n = 5, nburn_random = 100
real(kind=dp)      :: xx(n)
integer            :: ifix,offset
character (len=10) :: labels(0:1) = ["varying","fixed  "]
do offset=0,2
   do ifix=0,1
      call set_random_seed(fixed=(ifix == 1),offset=offset,nburn_random=nburn_random)
      call random_number(xx)
      write (*,"(a12,100f10.4)") trim(labels(ifix)),xx
   end do
end do
end program xset_random_seed