I run a fortran multithreaded code that uses ffts. I am now trying to switch to
fttw available thru mkl. When I insert the code below the results more
often than not are nans or infinity even though only the plan call is executed
and no actuall ffts are carried out using fftw. The standard suspect in such cases
would be out-of-bounds operations. I cannot at the moment see any out-of-bounds
problems when creating the plan. But here there is still the problem of the interface.
I wonder if i have set up things correcly. The compile and link calls are
ifort -O3 -r8 -openmp -fpp -parallel -mcmodel=medium -i-dynamic -shared-intel -mkl
This is run on a six core Linux machine.
Has anyone any suggestions as to what may be going on?
module FFTW3 use, intrinsic :: iso_c_binding include 'fftw3.f03' type(C_PTR), save :: plan_r2c real*8, allocatable, SAVE :: WWW1( :, : ), WWW2( :, : ) end module .... Use FFTW3 Integer :: Irank, ndim(1), InEmbeded(1), OutEmbeded(1) ... ! Note NX32 > N1_32 + 4 allocate( WWW1( NX32*M2, M3 ), STAT = IERR ) allocate( WWW2( NX32*M2, M3 ), STAT = IERR ) Irank = 1 ndim(1) = N1_32 InEmbeded(1) = NX32 OutEmbeded(1) = NX32 / 2 call dfftw_plan_many_dft_r2c( plan_r2c, Irank, ndim(1), M2*M3, WKX_1, InEmbeded(1), 1, & NX32, WKX_2, OutEmbeded(1), 1, NX32/2, FFTW_MEASURE )
The plan call is:
call dfftw_plan_many_dft_r2c( plan_r2c, Irank, ndim(1), M2*M3, WWW1, InEmbeded(1), 1, &
NX32, WWW2, OutEmbeded(1), 1, NX32/2, FFTW_MEASURE )
There is a test for the creation of the plan. See code extract below @ line 98.
I use fftw with a c code without any problem.
It is now that I am trying to insert fftw into a fortran code that all
hell broke loose.
I include the relevant subroutine (somewhat pruned) that contains the
relevant call. This is part of a much larger code. The problem is that just
this one call to the plan routine most times will contaminate with garbage other
parts of the code.
Any ideas would be appreciated.
module FFTW3 use, intrinsic :: iso_c_binding include 'fftw3.f03' type(C_PTR), save :: plan_r2c, plan_c2r real*8, allocatable, SAVE :: WWW1( :, : ), WWW2( :, : ) end module !=============================================== Subroutine Setup Use Tables Use FFTW3 Implicit None Integer :: I, J, K, JEL, KEL, Jump, II, JJ, KK Real*8 :: temp, VS_MAX, T1, T2 Integer :: Irank, ndim(1), InEmbeded(1), OutEmbeded(1) N1P1 = N1 + 1 N1P2 = N1 + 2 N1_32 = N1 + (N1P1/2) M2M1 = M2 - 1 M3M1 = M3 - 1 JLEN = N2*M2M1 - 1 KLEN = N3*M3M1 - 1 JFULL = N2*M2M1 + 1 KFULL = N3*M3M1 + 1 JPLEN = N2*M2M1 + 1 ! Pressure matrix y-size KPLEN = N3*M3M1 + 1 ! .. .. z-size RDT = 1.0D0 / DT RN1_32 = 1.0d0 / DFLOAT(N1_32) ! - Time staping factors EULER_1 = 1 EULER_2 = 2 NORDER = EULER_1 ! - Euler_1 settings. Alpha = RDT Beta = RDT Gamma = 0.0D0 Delta = 1.0D0 Epsilon = 0.0D0 TIME = 0.0D0 NTIME = 0 N_RANDOM = 0 T_RANDOM = 0 FlowRate_1 = 0.0D0 FlowRate_2 = 0.0D0 [ stuff delete ] T1 = pi*DFLOAT(N1)/XLEN Wv( N1+1 ) = T1 Wv_SQ( N1+1 ) = T1*T1 Call FFTFAX( N1, IFAX, TRIGS ) Call FFTFAX( N1_32, IFAX_32, TRIGS_32 ) !=========================================================== ! FFTW write(6,*)' ---- NX32, N1_32 = ', NX32, N1_32 allocate( WWW1( NX32*M2, M3 ), STAT = IERR ) if( Ierr .ne. 0 ) STOP ' ** FAILED WWWW1 ' allocate( WWW2( NX32*M2, M3 ), STAT = IERR ) if( Ierr .ne. 0 ) STOP ' ** FAILED WWWW2 ' Irank = 1 ndim(1) = N1_32 InEmbeded(1) = NX32 OutEmbeded(1) = NX32/2 call dfftw_plan_many_dft_r2c( plan_r2c, Irank, ndim(1), M2*M3, WWW1, InEmbeded(1), 1, & NX32, WWW2, OutEmbeded(1), 1, NX32/2, FFTW_MEASURE ) if( .not. c_associated(plan_r2c) ) STOP " *** FAILED TO CREAT FFTW PLAN" write(6,*) ' --- PLAN ', c_associated(plan_r2c) ! - Y,Z Grid and parameters. Call Grid Call ZWGLJD( YGL, RHY, M2, 0.0D0, 0.0D0 ) Call ZWGLJD( ZGL, RHZ, M3, 0.0D0, 0.0D0 ) Call DGLJD( DYGL, DYGLT, YGL, M2, MY, 0.0D0, 0.0D0 ) Call DGLJD( DZGL, DZGLT, ZGL, M3, MZ, 0.0D0, 0.0D0 ) ! - Set up filter data. Call ZWGLJD( YGL_F, WKP_1, M2_F, 0.0D0, 0.0D0 ) ! Y filter grid. WKP_1 is a dummy. Call ZWGLJD( ZGL_F, WKP_1, M3_F, 0.0D0, 0.0D0 ) ! Z filter grid. WKP_1 is a dummy. ! - Interpolation matrices : M2 -> M2_F and M3 -> M3_F Call IGLLMD( GLY2F, WKP_1, YGL, YGL_F, M2, M2_F, MY, MY, -999 ) Call IGLLMD( GLZ2F, WKP_1, ZGL, ZGL_F, M3, M3_F, MZ, MZ, -999 ) ! - Reverse interpolation : M2_F -> M2 and M3_F -> M3 Call IGLLMD( F2GLY, WKP_1, YGL_F, YGL, M2_F, M2, MY, MY, -999 ) Call IGLLMD( F2GLZ, WKP_1, ZGL_F, ZGL, M3_F, M3, MZ, MZ, -999 ) ! [stuff deleted] End Subroutine Setup
FFTW3 has Fortran interface. Why do you use C binding in your Fortran code? Please take a look at the MKL sample code that uses FFTW3 in Fortran code and see if it is helpful. The sample code is found in $MKLROOT/examples/fftw3xf/source/dp_plan_many_dft_r2c.f90
Thank you for your interest.
I have downloaded a number of test codes . All work ok.
My problem must be code specific.
I copied the module with the bindings from fthe fftw site. The c bindings
are meant to provide fortran types compatible with C which in this case
is the c_ptr. Anyway this is not the problem.
if i inderstood the mkl manual correctly, unlike the original fftw, the mkl
version is not thread safe even when different plans and arrays are used.
A single muntithreading call is ok tho.
This is a showstopper, unless u know different.
I don't know what problem you are talking about. If your question is about whether or not mkl functions are thread-safe then, just like Zhang has mentioned back in 2014, all mkl functions are thread-safe.