[Bug fortran/108336] New: Repeatable random_numbers with openmp

[Bug fortran/108336] New: Repeatable random_numbers with openmp

[gareth.davies.ga.code at gmail dot com]

https://gcc.gnu.org/bugzilla/show_bug.cgi?id=108336

            Bug ID: 108336
           Summary: Repeatable random_numbers with openmp
           Product: gcc
           Version: 11.3.0
            Status: UNCONFIRMED
          Severity: normal
          Priority: P3
         Component: fortran
          Assignee: unassigned at gcc dot gnu.org
          Reporter: gareth.davies.ga.code at gmail dot com
  Target Milestone: ---

The aim is to get repeatable random numbers with openmp. This report was
derived from a discussion including input from gfortran developer Steve Kargl:
https://fortran-lang.discourse.group/t/replicable-monte-carlo-simulation-with-openmp/4972/19

AFAIK the following program should give repeatable random numbers using openmp.
But it doesn't - there are occasional changes in results, suggestive of a race
condition. 


$ gfortran --version
GNU Fortran (Ubuntu 11.3.0-1ubuntu1~22.04) 11.3.0
Copyright (C) 2021 Free Software Foundation, Inc.
This is free software; see the source for copying conditions.  There is NO
warranty; not even for MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.

$ cat random_number7.f90 
program main
   use omp_lib
   implicit none
   integer :: n1, i, n, nt
   real :: x1, x2
   real, allocatable :: threadval(:)

   n1 = 2

   ! Initialize the master state to a repeatable sequence.
   call random_init(.true., .false.)

   ! Store final random number on each thread
!$omp parallel
   nt = omp_get_num_threads()
!$omp end parallel
   allocate(threadval(nt))
   threadval = -1.0

!$omp parallel private(x1, x2, i), shared(n1, threadval)

   ! Trigger non-master threads to get their unique seed
   if(omp_get_thread_num() /= 0) call random_number(x1)
!$omp barrier
   ! Get some random numbers -- almost repeatable (but not)
!$omp do
   do i = 1, n1
      call random_number(x1)
      call random_number(x2)
   end do
!$omp end do
   ! Store final RN on each thread
   threadval(omp_get_thread_num()+1) = x2 

!$omp end parallel

   print*, minval(threadval), maxval(threadval)
end program

$ gfortran -fopenmp random_number7.f90 -o random_number7
$ for i in {1..10}; do OMP_NUM_THREADS=6 ./random_number7; done
  0.191325366      0.876643896    
  0.191325366      0.480524838    
  0.191325366      0.876643896    
  0.191325366      0.804360509    
  0.191325366      0.804360509    
  0.191325366      0.480524838    
  0.191325366      0.670394361    
  0.191325366      0.804360509    
  0.157898605      0.191325366    
  0.191325366      0.670394361

[Bug fortran/108336] Repeatable random_numbers with openmp

[pinskia at gcc dot gnu.org]

https://gcc.gnu.org/bugzilla/show_bug.cgi?id=108336

--- Comment #1 from Andrew Pinski <pinskia at gcc dot gnu.org> ---
Most likely there is only one seed per program so the random # generator is
using locks. So the order is dependent on which thread calls the function.

Now having a seed per thread might be a good idea but it does have other kind
of issues. Even initializing it and even Making sure it does not reinitialized
each time a new openmp thread is launched.

[Bug fortran/108336] Repeatable random_numbers with openmp

[gareth.davies.ga.code at gmail dot com]

https://gcc.gnu.org/bugzilla/show_bug.cgi?id=108336

--- Comment #2 from Gareth Davies <gareth.davies.ga.code at gmail dot com> ---
(In reply to Andrew Pinski from comment #1)
> Most likely there is only one seed per program so the random # generator is
> using locks. So the order is dependent on which thread calls the function.
> 
> Now having a seed per thread might be a good idea but it does have other
> kind of issues. Even initializing it and even Making sure it does not
> reinitialized each time a new openmp thread is launched.

There is some documentation as to how it's supposed to work here:
https://gcc.gnu.org/onlinedocs/gfortran/RANDOM_005fSEED.html

It sounds like we should be able to get repeated random numbers - even if their
mapping to the threads does vary. Indeed, the results are pretty close to that.

[Bug libfortran/108336] Repeatable random_numbers with openmp

[pinskia at gcc dot gnu.org]

https://gcc.gnu.org/bugzilla/show_bug.cgi?id=108336

--- Comment #3 from Andrew Pinski <pinskia at gcc dot gnu.org> ---
Yes the race condition is the following:
the order of calling random_number on which thread first.

init_rand_state definitely takes a lock when reading/writing from/to
master_state.

I see no other race condition in the code.
The code is in intrinsics/random.c 

https://gcc.gnu.org/git/?p=gcc.git;a=blob;f=libgfortran/intrinsics/random.c;h=b5732e6c25a14359745c5f698c2dd4cab6f3c81f;hb=HEAD

https://gcc.gnu.org/git/?p=gcc.git;a=blob;f=libgfortran/intrinsics/random_init.f90;h=c2a6a0db9ee34c968494d1d2692ac4a48097251f;hb=HEAD

[Bug libfortran/108336] Repeatable random_numbers with openmp

[anlauf at gcc dot gnu.org]

https://gcc.gnu.org/bugzilla/show_bug.cgi?id=108336

--- Comment #4 from anlauf at gcc dot gnu.org ---
There is one flaw with the testcase: when OMP_NUM_THREADS > n1,
array elements threadval(n1+1:OMP_NUM_THREADS) are filled with
undefined values.

When I replace the line

   if(omp_get_thread_num() /= 0) call random_number(x1)

by the following loop:

   do i = 1, nt
      if (i > 1 .and. omp_get_thread_num() == i-1) call random_number(x1)
!$omp barrier
   end do

and furthermore do

   print*, minval(threadval(1:n1)), maxval(threadval(1:n1))

I get repeatable results.

Am I missing something?

[Bug libfortran/108336] Repeatable random_numbers with openmp

[pinskia at gcc dot gnu.org]

https://gcc.gnu.org/bugzilla/show_bug.cgi?id=108336

--- Comment #5 from Andrew Pinski <pinskia at gcc dot gnu.org> ---
Maybe instead of:

!$omp parallel
   nt = omp_get_num_threads()
!$omp end parallel

Do:
  nt = omp_get_max_threads()

that should improve the situtation of not allocating enough for the array.

[Bug libfortran/108336] Repeatable random_numbers with openmp

[gareth.davies.ga.code at gmail dot com]

https://gcc.gnu.org/bugzilla/show_bug.cgi?id=108336

--- Comment #6 from Gareth Davies <gareth.davies.ga.code at gmail dot com> ---
(In reply to anlauf from comment #4)
> There is one flaw with the testcase: when OMP_NUM_THREADS > n1,
> array elements threadval(n1+1:OMP_NUM_THREADS) are filled with
> undefined values.

Yes, my bad in trying to simplify the test case...

We should have n1 >= nt.

> 
> When I replace the line
> 
>    if(omp_get_thread_num() /= 0) call random_number(x1)
> 
> by the following loop:
> 
>    do i = 1, nt
>       if (i > 1 .and. omp_get_thread_num() == i-1) call random_number(x1)
> !$omp barrier
>    end do
> 
> and furthermore do
> 
>    print*, minval(threadval(1:n1)), maxval(threadval(1:n1))
> 
> I get repeatable results.
> 
> Am I missing something?

Thanks, I also get repeatable results with those changes. 

A "weaker" approach to getting deterministic results without forcing a 1:1
mapping of random numbers to threads is to replace:

>    do i = 1, nt
>       if (i > 1 .and. omp_get_thread_num() == i-1) call random_number(x1)
> !$omp barrier
>    end do

with:

!$omp critical
   if(omp_get_thread_num() /= 0 ) call random_number(x1)
!$omp end critical 
!$omp barrier

This is what I was looking for.