Replicable Monte Carlo simulation with OpenMP

Help
3

Another issue in your code is the race condition on x1 and x2: they have to be private (as in the code proposed by @gareth )

I would initialize the seeds with deterministic random numbers as well…

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

  n1 = 1000*1000

  call random_seed0(size=n)
  seed0 = [ (i, i=0, n-1) ]
  call random_seed(put=seed0)

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

  !$omp single
  allocate( fseed(n,omp_get_num_threads()) )
  call random_number(fseed)
  seed = floor(fseed*huge(seed))
  !$omp end single

  ! Deterministic seed, distinct on each thread.
  call random_seed(put=seed(:,omp_get_thread_num()))

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

That said, are random_seed() and random_number() garanteed to be thread-safe ?