I have trouble using the fftw_plan routines for 2D datas in an openMP environment. My compiler is gfortran, and my OS is Ubuntu 12.04. I compiled fftw with the --enable-openmp option, and the library seems to be linked correctly. However, when I get into trouble when I run the following minimal exemple:
PROGRAM test_fft_OMP
USE, intrinsic :: iso_c_binding
USE OMP_LIB
IMPLICIT NONE
INCLUDE 'fftw3.f03'
INTEGER*4, PARAMETER :: N=512
REAL(C_DOUBLE), DIMENSION(N) :: dbgr
COMPLEX(C_DOUBLE_COMPLEX), DIMENSION(N/2+1) :: dbgc
COMPLEX(C_DOUBLE_COMPLEX), DIMENSION(N) :: dbgcc
COMPLEX(C_DOUBLE_COMPLEX), DIMENSION(N,N) :: D2CA
COMPLEX(C_DOUBLE_COMPLEX), DIMENSION(N,N) :: D2CB
REAL(C_DOUBLE), DIMENSION(N,N) :: D2R
COMPLEX(C_DOUBLE_COMPLEX), DIMENSION(N/2+1,N):: D2C
INTEGER :: nthreads, void
TYPE(C_PTR) :: plan
INTEGER(C_INT), DIMENSION(2) :: n2 = [N,N]
INTEGER(C_INT), DIMENSION(1) :: n1 = [N]
nthreads = omp_get_num_threads()
void = fftw_init_threads()
CALL fftw_plan_with_nthreads(nthreads)
plan=fftw_plan_dft_c2r_1d(N, dbgc , dbgr,FFTW_MEASURE)
!plan=fftw_plan_dft_c2r_2d(N,N, D2C , D2R ,FFTW_MEASURE)
plan=fftw_plan_dft_1d(N, dbgcc , dbgcc,FFTW_FORWARD,FFTW_MEASURE)
plan=fftw_plan_dft_2d(N,N, D2CA , D2CA ,FFTW_FORWARD,FFTW_MEASURE)
plan=fftw_plan_dft(1,n1, dbgc , dbgc ,FFTW_FORWARD,FFTW_MEASURE)
plan=fftw_plan_dft(2,n2, D2CA , D2CA ,FFTW_FORWARD,FFTW_MEASURE)
!plan=fftw_plan_dft_2d(N,N, D2CA , D2CB ,FFTW_FORWARD,FFTW_MEASURE)
!plan=fftw_plan_dft(2,n2, D2CA , D2CB ,FFTW_FORWARD,FFTW_MEASURE)
END PROGRAM test_fft_OMP
The code compilation is fine, using
gfortran -fopenmp -o execmp test.f90 -I/usr/local/lib/ -L/usr/local/lib/ -lfftw3_omp -lfftw3 -lm -lpthread
However, if I uncomment any of the three !plan lines, I obtain a segmentation fault if I run the executable.
If I summarize my problem:
I can create plans for 1D transforms, or 2D in place transforms
I get a seg fault for 2D out of place transforms, whatever the syntax used
The FFTW planner crashes for large arrays (hence my problem with 2D arrays, but only above ~512x512) if the arrays are declared statically instead of using ALLOCATE for a dynamic memory allocation.
More specifically, this variable declaration makes the program crash when fftw_plan is invoked,
COMPLEX(C_DOUBLE_COMPLEX), DIMENSION(N,N) :: D2C
whereas this one works fine:
COMPLEX(C_DOUBLE_COMPLEX), ALLOCATABLE, DIMENSION(:,:) :: D2C
ALLOCATE(D2C(N,N))