Initializing RANDOM_NUMBER

1

For codes using RANDOM_NUMBER, you may wish to call RANDOM_SEED to have different variates generated on each program run or to get the same variates, for reproducibility. I came up with this:

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
      write (*,*) "in set_random_seed, nseeds, nseeds_max =",nseeds,nseeds_max,"need nseeds <= nseeds_max"
      stop
   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 = 0
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

Running with gfortran, I get

c:\fortran\2015>a.exe
     varying    0.6095    0.0214    0.8085    0.5975    0.0457
       fixed    0.6144    0.9234    0.1654    0.1935    0.0150
     varying    0.8507    0.1854    0.3608    0.7726    0.6638
       fixed    0.3029    0.6460    0.0766    0.8636    0.6028
     varying    0.9658    0.2016    0.5814    0.6148    0.9721
       fixed    0.9545    0.7794    0.6587    0.9352    0.1320

c:\fortran\2015>a.exe
     varying    0.5870    0.2881    0.8306    0.7200    0.3491
       fixed    0.6144    0.9234    0.1654    0.1935    0.0150
     varying    0.6186    0.9283    0.0240    0.0694    0.4192
       fixed    0.3029    0.6460    0.0766    0.8636    0.6028
     varying    0.3099    0.7066    0.9274    0.7142    0.8246
       fixed    0.9545    0.7794    0.6587    0.9352    0.1320

c:\fortran\2015>a.exe
     varying    0.4781    0.9395    0.5838    0.5719    0.4711
       fixed    0.6144    0.9234    0.1654    0.1935    0.0150
     varying    0.8101    0.2934    0.5815    0.6075    0.2483
       fixed    0.3029    0.6460    0.0766    0.8636    0.6028
     varying    0.5737    0.8283    0.0903    0.5931    0.1938
       fixed    0.9545    0.7794    0.6587    0.9352    0.1320

c:\fortran\2015>a.exe
     varying    0.2166    0.5586    0.3375    0.3679    0.8439
       fixed    0.6144    0.9234    0.1654    0.1935    0.0150
     varying    0.8032    0.6942    0.0409    0.8707    0.5342
       fixed    0.3029    0.6460    0.0766    0.8636    0.6028
     varying    0.1004    0.9886    0.4712    0.6127    0.0289
       fixed    0.9545    0.7794    0.6587    0.9352    0.1320

so I get either the same or different variates, as desired. I get similar results with Intel Fortran. Initially, I did not burn off variates (corresponding to nburn_random=0 in the code above), but then the random streams start out the same for small changes in offsets and take some time to diverge, for example

 varying    0.9703    0.6780    0.9745    0.7573    0.6670
   fixed    0.4593    0.8546    0.7349    0.9362    0.0186
 varying    0.9703    0.6781    0.4471    0.6909    0.8093
   fixed    0.4593    0.8546    0.9107    0.8151    0.9856
 varying    0.9703    0.6781    0.6229    0.4272    0.6665
   fixed    0.4593    0.8546    0.0865    0.9029    0.1284

Gfortran, Intel, and g95 use seed sizes of 8, 2, and 4. Does any compiler using a seed size larger than 8? Suggested improvements to the code are welcome.

1 Like
2

The NAG Fortran compiler needs 630.

And you need to be careful with your "seed = " statements, they reallocate to the size of the RHS.

1 Like
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
4

seed array reallocation on assignment still produces arrays too small for PUT=.

6

gfortran older than some version, which you can look up via release notes, uses 12.

Thanks – now my program below hard codes 12 random integers. I hope I also fixed the bug reported by @themos. Could someone with a NAG compiler test it? I never intentionally use allocation upon assignment to change the size of an allocatable array, because of the kind of bugs I have produced in this thread.

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 = 12, seeds_vec(nseeds_max) = [824514020,218904035,384790913, &
                                 510442021,939900036,939295695,826403327,935378387, &
                                 634734772,413176190,91069182,818551196]
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
   else
      seed = seeds_vec(:nseeds) + offset_
   end if
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
7 
   seed(1:nseeds_max) = xseed(nburn+1:) * 1000000000 + offset_
   seed(nseeds_max+1:) = 0
9

size(RHS)=12, size(LHS)=630, seed was resized to 12.

10

I hope this is right – tested with gfortran and Intel Fortran.

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 = 12, seeds_vec(nseeds_max) = [824514020,218904035,384790913, &
                                 510442021,939900036,939295695,826403327,935378387, &
                                 634734772,413176190,91069182,818551196]
integer                       :: i,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
   else
      seed = seeds_vec(:nseeds) + offset_
   end if
else ! set seed based on current time
   call system_clock(itime)
   seed = itime
   call random_seed(put=seed)
   call random_number(xseed)
   do i=1,nseeds
      if (i <= nseeds_max) then
         seed(i) = xseed(nburn+i)*1000000000 + offset_
      else
         seed(i) = 0
      end if
   end do
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
1 Like
11

Does the new RANDOM_INIT intrinsic help with this issue?

12

I just tried it. For a program lightly modified from the Fortran Wiki

program test_random_seed
implicit none
integer, parameter :: n = 5
real x(n), y(n)
call random_init(repeatable=.true.,image_distinct=.true.)
call random_number(x)
call random_init(repeatable=.true.,image_distinct=.true.)
call random_number(y)
if (any(x /= y)) stop "x(:) /= y(:)"
write (*,"(100f9.6)") x
end program test_random_seed

gfortran prints different values of x on each run, so the random numbers are repeatable within a program but not across runs. Intel Fortran build 20201112_000000 gives output

0.000000 0.025480 0.352516 0.666914 0.963055

each time. Can I use RANDOM_INIT to get the same random numbers across runs for both compilers? To get different random numbers for each run?

It also concerns me that for Intel the first random number produced is zero, although that is standard-conforming.