Replicable Monte Carlo simulation with OpenMP

Help
2

See documentation for random_init here. It doesn’t do what you expect. For example, .... The sequence of random numbers is different for repeated execution of the program.

Another factor is that random_init is designed for coarray parallelism - which is not the same as openmp parallelism.

The program below shows one way to get repeatable results, but only if OMP_NUM_THREADS is fixed. However the approach to setting the seed might be unwise.

program main
  use omp_lib
  implicit none
  integer :: n1, n2, i, n
  integer, allocatable :: seed(:)
  real(8) :: pi, x1, x2

  n1 = 1000*1000

  ! Make space for a deterministic seed
  call random_seed(size=n)
  allocate(seed(n))

  !$omp parallel firstprivate(seed), private(x1, x2, i), shared(n1), reduction(+:n2)
  n2 = 0

  ! Deterministic seed, distinct on each thread.
  seed = [ (i, i=0, size(seed)-1) ] + omp_get_thread_num()
  call random_seed(put=seed)

  !$omp do
  do i = 1, n1
    call random_number(x1)
    call random_number(x2)
    if (x1**2 + x2**2 <= 1.0_8) n2 = n2 + 1
  end do
  !$omp end do

  !$omp end parallel
  pi = 4.0_8*dble(n2)/dble(n1)
  print *, pi
end program

While this seems to work, I don’t know if we can be confident that these seeds will produce effectively independent random numbers. That probably depends on the random number generator.

1 Like