fortran+fftw=nan (sometimes)

[a_b_1]

Hello,

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?

Thank you.

--

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 )

[a_b_1]

Ooops CORRECTION.

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 )

[Zhang_Z_Intel]

After the plan call returns, did you check if plan_r2c is zero? Probably, the plan was not created successfully for some reason. Can you attach your full test code?
 

[a_b_1]

Thank you.

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

[Zhang_Z_Intel]

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

[a_b_1]

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.

--

[Zhang_Z_Intel]

All MKL functions, including FFT, are thread safe. It is clearly stated in the MKL User's Guide: http://software.intel.com/en-us/mkl_11.1_ug_lin

[Osinsky__Alexander]

I have literally the same problem. Has it been solved since then?

[Khang_N_Intel]

Hi Alexander,

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.

Best Regards,

Khang