使用Fortran的OpenMP FFTW不是线程安全的
我尝试在openMP和Fortran中使用快速傅立叶变换，但在并行执行时得到错误的结果，这也会在每个执行步骤中改变它们的值，在并行出错时显示典型的行为。
我的目标是一个简单的3d实数到复杂的转换。按照FFTW教程，除了对dfftw_execute_dft_r2c()的调用之外，我把所有的东西都移出了并行区域，但它似乎不起作用。
我使用配置了./configure --enable-threads --enable-openmp --enable-mpi的FFTW3.3.8，并使用gfortran program.f03 -o program.o -I/usr/include -fopenmp -lfftw3_omp -lfftw3 -g -Wall编译我的代码。
下面是我的程序的样子：
program use_fftw
use,intrinsic :: iso_c_binding
use omp_lib
implicit none
include 'fftw3.f03'
integer, parameter :: dp=kind(1.0d0)
integer, parameter :: Nx = 10
integer, parameter :: Ny = 5
integer, parameter :: Nz = 5
real(dp), parameter :: pi = 3.1415926d0
real(dp), parameter :: physical_length_x = 20.d0
real(dp), parameter :: physical_length_y = 10.d0
real(dp), parameter :: physical_length_z = 10.d0
real(dp), parameter :: lambda1 = 0.5d0
real(dp), parameter :: lambda2 = 0.7d0
real(dp), parameter :: lambda3 = 0.9d0
real(dp), parameter :: dx = physical_length_x/real(Nx,dp)
real(dp), parameter :: dy = physical_length_y/real(Ny,dp)
real(dp), parameter :: dz = physical_length_z/real(Nz,dp)
integer :: void, nthreads
integer :: i, j, k
real(dp):: d
complex(dp), allocatable, dimension(:,:,:) :: arr_out
real(dp), allocatable, dimension(:,:,:) :: arr_in
integer*8 :: plan_forward
allocate(arr_in( 1:Nx, 1:Ny, 1:Nz)); arr_in = 0
allocate(arr_out(1:Nx/2+1, 1:Ny, 1:Nz)); arr_out = 0
!------------------------------
! Initialize fftw stuff
!------------------------------
! Call before any FFTW routine is called outside of parallel region
void = fftw_init_threads()
if (void==0) then
write(*,*) "Error in fftw_init_threads, quitting"
stop
endif
nthreads = omp_get_num_threads()
call fftw_plan_with_nthreads(nthreads)
! plan execution is thread-safe, but plan creation and destruction are not:
! you should create/destroy plans only from a single thread
call dfftw_plan_dft_r2c_3d(plan_forward, Nx, Ny, Nz, arr_in, arr_out, FFTW_ESTIMATE)
!--------------------------------
! Start parallel region
!--------------------------------
!$OMP PARALLEL PRIVATE( j, k, d)
! Fill array with wave
! NOTE: wave only depends on x so you can plot it later.
!$OMP DO
do i = 1, Nx
d = 2.0*pi*i*dx
do j = 1, Ny
do k = 1, Nz
arr_in(i,j,k) = cos(d/lambda1)+sin(d/lambda2)+cos(d/lambda3)
enddo
enddo
enddo
!$OMP END DO
call dfftw_execute_dft_r2c(plan_forward, arr_in, arr_out)
!$OMP END PARALLEL
!-----------------
! print results
!-----------------
do i=1, Nx/2+1
do j=1, Ny
do k=1, Nz
write(*,'(F12.6,A3,F12.6,A3)',advance='no') real(arr_out(i,j,k)), " , ", aimag(arr_out(i,j,k)), " ||"
enddo
write(*,*)
enddo
write(*,*)
enddo
deallocate(arr_in, arr_out)
! destroy plans is not thread-safe; do only with single
call dfftw_destroy_plan(plan_forward)
end program use_fftw我也尝试了将FFTW (void = fftw_init_threads(); call fftw_plan_with_nthreads(nthreads); call dfftw_plan_dft_r2c_3d(...) )的初始化部分移动到并行区，使用!$OMP SINGLE块并随后与屏障同步，但情况并没有改善。
有谁可以帮我？
编辑：我可以在另一个系统上测试我的程序，问题仍然存在。所以问题显然不在我的openmp或FFTW的实现中，而是在程序本身的某个地方。
回答 1
Stack Overflow用户
发布于 2018-06-09 04:25:27
通常应在parallel区域的外部调用fftw执行例程。它们在它们内部有自己的并行区域，它们将负责在规划期间与您请求的那么多线程并行运行转换。它们将重用您现有的OpenMP线程。
你也可以在一个并行的区域内调用它们，但是在不同的数组上，而不是在相同的数组上！然后你的计划应该计划使用一个线程。然后，例如，每个线程将执行数组切片的2D变换。
线程安全意味着您可以并发地调用这些例程，但每个例程针对不同的数据。
https://stackoverflow.com/questions/50765932
复制相似问题