Turn on suggestions

Auto-suggest helps you quickly narrow down your search results by suggesting possible matches as you type.

Showing results for

- Intel Community
- Software
- Software Development SDKs and Libraries
- Intel® oneAPI Math Kernel Library
- fortran+fftw=nan (sometimes)

- Subscribe to RSS Feed
- Mark Topic as New
- Mark Topic as Read
- Float this Topic for Current User
- Bookmark
- Subscribe
- Mute
- Printer Friendly Page

a_b_1

Beginner

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

04-02-2014
11:51 AM

223 Views

fortran+fftw=nan (sometimes)

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 )

Link Copied

8 Replies

a_b_1

Beginner

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

04-02-2014
11:58 AM

223 Views

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

Employee

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

04-02-2014
04:17 PM

223 Views

a_b_1

Beginner

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

04-03-2014
12:51 PM

223 Views

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

Employee

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

04-04-2014
04:03 PM

223 Views

a_b_1

Beginner

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

04-07-2014
01:56 AM

223 Views

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

Employee

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

04-07-2014
10:30 AM

223 Views

Osinsky__Alexander

Beginner

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

11-06-2018
08:04 AM

223 Views

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

Khang_N_Intel

Employee

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

11-07-2018
10:48 AM

223 Views

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

Topic Options

- Subscribe to RSS Feed
- Mark Topic as New
- Mark Topic as Read
- Float this Topic for Current User
- Bookmark
- Subscribe
- Printer Friendly Page

For more complete information about compiler optimizations, see our Optimization Notice.