Thread safe FFTW

Hi,

I am attempting to incorporate a layer of OpenMP code within an MPI program.
Unfortunately I seem to be falling foul of thread safety when attempting to perform FFTs within the OpenMP threads.

I am fairly convinced that thread safety is the issue as the code executes without issue if the call to dfftw_execute_dft is performed within an omp critical region but fails in parallel with a double free or corruption error.

As far as I can tell the code is structured in the manner advised by the fftw documentation: The plan is generated before the OpenMP region is created and dfftw_execute_dft is supposed to be a thread safe routine for executing the FFTs.

I have attempted to use the FFTW_UNALIGNED flag as specified in a post that describes a very similar problem on: http://fortran.itgroups.info/69/9/thread-1901648.html with no joy.

Before I resort to changing the manner in which the workload is distributed I wondered of there may be anything I missed?

Thanks in advance,

Karl

Hi Karl,

I’m not sure. Using the sample code in the link you reference, if I move the create and destroy plan calls outside the OpenMP region, the code runs correctly. I did make the mistake of forgetting to move the “plan_r2c” out of the private clause which caused a segv, but was find once made shared. Could this be your issue as well?

  • Mat
% cat main.F 
program main
implicit none
#include "fftw3.inc"
integer nx,ny
parameter(nx=128,ny=8)
double precision val(nx,ny)
double precision tp(nx), tf(nx)
integer i,j,k
real*8 deltax,rlenx
real*8 pi,amp

integer*8 plan_r2c
integer tid,OMP_GET_THREAD_NUM

pi=4.d0*atan(1.d0)
rlenx=2.d0*pi
deltax=rlenx/dble(nx)
! write(*,*)'rlenx:',rlenx

amp=2.d0

do j=1,ny
do i=1,nx
val(i,j)=amp*dsin(2.d0*pi*2.d0*deltax*dble(i-1)/rlenx)
enddo
enddo
call dfftw_plan_dft_r2c_1d(plan_r2c,nx,tp,tf,FFTW_MEASURE)

!CALL OMP_SET_NUM_THREADS(4)

!omp parallel do private(tid,k,plan_r2c,tp,tf)
!$omp parallel do private(tid,k,tp,tf)
!omp parallel do default(private)
do k =1,ny
tid = OMP_GET_THREAD_NUM()
write(*,*) 'Thread ID = ', tid, ' Loop # = ', k

do i=1,nx
tp(i)=val(i,k)
tf(i)=dble(0.0)
enddo

!call dfftw_plan_dft_r2c_1d(plan_r2c,nx,tp,tf,FFTW_MEASURE)
call dfftw_execute(plan_r2c)
!call dfftw_destroy_plan(plan_r2c)

write(*,*) 'plan_r2c = ', plan_r2c

enddo
call dfftw_destroy_plan(plan_r2c)

end
% pgf90 main.F -Mfree -L/local/home/colgrove/fftw/lib -lfftw3 -g -mp
% setenv OMP_NUM_THREADS 4
% a.out
 Thread ID =             0  Loop # =             1
 Thread ID =             2  Loop # =             5
 plan_r2c =                   8760496
 Thread ID =             0  Loop # =             2
 plan_r2c =                   8760496
 Thread ID =             2  Loop # =             6
 Thread ID =             3  Loop # =             7
 plan_r2c =                   8760496
 plan_r2c =                   8760496
 plan_r2c =                   8760496
 Thread ID =             3  Loop # =             8
 plan_r2c =                   8760496
 Thread ID =             1  Loop # =             3
 plan_r2c =                   8760496
 Thread ID =             1  Loop # =             4
 plan_r2c =                   8760496