Fortran结合MPI实现有放回重采样的独立随机数问题

使用Fortran搭配MPI对规模为$10^7$的数组DATA执行有放回重采样，需完成1000次重采样并分配到10个处理器上。当前代码中所有处理器调用random_number生成的new_index_real数组完全相同，无法保证1000次重采样彼此独立。相关代码片段如下：

real :: DATA(10^7), DATA_resample(10^7)
integer :: org_index_int(10^7)
integer :: new_index_int(10^7)
real :: new_index_real(10^7)

call MPI_INIT(ierror)
call MPI_COMM_RANK(MPI_COMM_WORLD, mid, ierror)

call random_init(.true., .true.)

! resample with replacement
! for myid=0 : irep_beg, irep_end = 1, 100
! for myid=9 : irep_beg, irep_end = 901, 1000
do irep=irep_beg, irep_end
   call random_number(new_index_real)
   new_index_int=floor(new_index_real*10^7)+1
   DATA_resample=DATA(new_index_int)
end do

call MPI_FINALIZE(ierror)

解决方案 #

1. 为每个MPI进程设置独立的随机数种子 #

所有进程默认使用相同的随机数种子初始化，导致random_number输出完全一致。需利用进程编号（mid）为每个进程生成差异化种子：

integer :: seed_size, i
integer, allocatable :: seed(:)

call MPI_INIT(ierror)
call MPI_COMM_RANK(MPI_COMM_WORLD, mid, ierror)

! 手动初始化随机数种子，基于进程号差异化
call random_seed(size=seed_size)
allocate(seed(seed_size))
call random_seed(get=seed)
! 给每个进程的种子添加进程号相关偏移，确保种子不重复
seed = seed + mid * 100000
call random_seed(put=seed)
deallocate(seed)

若使用支持Fortran 2008及以上的编译器，可简化为：

call random_init(.true., .true.)  ! 第二个参数设为.true.时，不同MPI进程会获得独立种子

2. 修复数组大小的语法错误 #

Fortran中^是按位异或运算符，而非幂运算，10^7会被解析为13，完全不符合数组规模需求。需改为10**7表示10的7次方：

real :: DATA(10**7), DATA_resample(10**7)
integer :: org_index_int(10**7)
integer :: new_index_int(10**7)
real :: new_index_real(10**7)

同时注意new_index_int=floor(new_index_real*10^7)+1中的10^7也需改为10**7，否则生成的索引范围错误。

3. 进程内重采样的独立性保障 #

每个进程内的do循环中，每次调用random_number会自动推进随机数生成器的内部状态，因此循环内的每一次重采样都会得到独立的随机索引，无需额外操作。

内容的提问来源于stack exchange，提问作者Jia