在Fortran中，有没有一种用随机数生成器进行并发计算的简单方法？

Question

我有一些我用Fortran 95编写的有限元代码，我已经对它们进行了优化，现在我可以得到16Mil以上的数据。在2 2GB内存占用下工作的元素。

我的代码的源函数不是平滑的，所以我使用(分层的)蒙特卡洛方法进行积分，这需要一个随机数生成器来选择样本位置

我尝试过使用-fopenmp -Ofast -ftree-parallelize-loops=4使用gfortran-9进行编译，但是使用随机数生成器的循环不能并行运行。我尝试过do concurrent，但显然不起作用，因为random_number不是“纯粹的”。https://stackoverflow.com/a/32637737/2372254我也试过阻塞我的循环，但也不起作用。

下面是我正在讨论的代码

        do  k=1,n_els ! total elements is n_els**2. This is block
            do i=1+ (k-1)*n_els ,k*n_els
                supp_vec = 0 
                integ_vec = 0.0_wp

                ! in this subroutine I call random_number
                call do_element(ind, n_els, i, num_points_per_strat, &
                                strat_rows, strat_cols, supp_vec, integ_vec) 
                do j=1, 4
                    sc_vec(supp_vec(j) ) = integ_vec(j)
                end do
                ! give some info about progress

                if (mod( i , (n_els**2)/10) == 0) print*, i*10/((n_els**2)/10), "% done"
            end do
        end do

似乎我可以将块写到一个文件中，然后调用n例程的不同实例。我想一定有一种更干净的方法来做到这一点。有什么建议可以让它更快地进行吗？

我正在考虑首先将一个块大小的点(取决于内存限制)写到一个数组中，然后在子例程调用中提供它。在我尝试之前，我想我会看看是否有人有更好的方法的建议。在可能的情况下尽量减少内存占用将是一件很好的事情。

Answer 1

在版本7和更新版本中，GFortran有一个并行的随机数生成器。在实现它时，下面是我用来验证性能是否确实随着线程数量的增加而扩展的OpenMP代码(来自https://gcc.gnu.org/ml/gcc-patches/2015-12/msg02110.html )：

! Benchmark generating random numbers
! Janne Blomqvist 2015
program randbench
#ifdef _OPENMP
  use omp_lib
#endif
  implicit none
  integer, parameter :: dp=kind(0.d0)          ! double precision
  integer, parameter :: i64 = selected_int_kind(18) ! At least 64-bit integer

#ifdef _OPENMP
  print *, "Using up to ", omp_get_max_threads(), " threads."
#endif
  call genr4
  call genr8
contains
  subroutine genr4
    integer, parameter :: n = int(1e7)
    real, save :: r(n)
    integer :: i
    integer(i64) :: t1, t2, td
#ifdef _OPENMP
    integer :: blocks, blocksize, l, h
#endif
    Print *, "Generate default real random variables"
    call system_clock (t1)
    !$omp parallel do private(i)
    do i = 1, n
       call random_number(r(i))
    end do
    !$omp end parallel do
    call system_clock (t2)
    td = t2 - t1
    print *, "Generating ", n, " default reals individually took ", td, " ticks."
    call system_clock (t1)
#ifdef _OPENMP
    blocks = omp_get_max_threads()
    blocksize = n / blocks
    !$omp parallel do private(l,h,i)
    do i = 0, blocks - 1
       l = i * blocksize + 1
       h = l + blocksize - 1
       !print *, "Low: ", l, " High: ", h
       call random_number(r(l:h))
    end do
#else
    call random_number(r)
#endif
    Call system_clock (t2)
    print *, "Generating ", n, " default reals as an array took  ", t2-t1, &
         " ticks. => ind/arr = ", real(td, dp) / (t2-t1)
  end subroutine genr4

  subroutine genr8
    integer, parameter :: n = int(1e7)
    real(dp), save :: r(n)
    integer :: i
    integer(i64) :: t1, t2, td
#ifdef _OPENMP
    integer :: blocks, blocksize, l, h
#endif
    print *, "Generate double real random variables"
    call system_clock (t1)
    !$omp parallel do
    do i = 1, n
       call random_number(r(i))
    end do
    call system_clock (t2)
    td = t2 - t1
    print *, "Generating ", n, " double reals individually took  ", td, " ticks."
    call system_clock (t1)
#ifdef _OPENMP
    blocks = omp_get_max_threads()
    blocksize = n / blocks
    !$omp parallel do private(l,h,i)
    do i = 0, blocks - 1
       l = i * blocksize + 1
       h = l + blocksize - 1
       !print *, "Low: ", l, " High: ", h
       call random_number(r(l:h))
    end do
#else
    call random_number(r)
#endif
    call system_clock (t2)
    print *, "Generating ", n, " double reals as an array took   ", t2-t1, &
         " ticks. => ind/arr = ", real(td, dp) / (t2 -t1)
  end subroutine genr8

end program