Intel® MPI Library
Get help with building, analyzing, optimizing, and scaling high-performance computing (HPC) applications.

3-D FFT via 1-D conversion

Masrul
Beginner
662 Views

Hello all,

I am trying to do 3-D FFT complex to complex. First i converted Q to complex QC then reshape 3D to 1D. Then I did FFT. Could anybody tell, if i did any mistake?

            ! convert Q to complex QC
            QC=cmplx(Q,0.d0,16)

           ! Reshape 3D to 1D
            QC_1d=reshape(QC,(/nfftdim1*nfftdim2*nfftdim3/))

            ! Initializing  DFFT stuff
            FFT_length(1)=nfftdim1
            FFT_length(2)=nfftdim2
            FFT_length(3)=nfftdim3
          

            ! Initializing  Forward DFFT
            status=DftiCreateDescriptor(forward_desc,DFTI_DOUBLE,DFTI_COMPLEX, 3, FFT_length)
            status=DftiCommitDescriptor(forward_desc)


            ! Doing Forward FFT
            status=DftiComputeForward(forward_desc, QC_1d)

 

0 Kudos
6 Replies
McCalpinJohn
Honored Contributor III
662 Views

A three-dimensional FFT is not the same as a one-dimensional FFT on the concatenated data.

It is common to decompose large one-dimensional FFTs into pseudo-2D or pseudo-3D arrays, but the transforms have to be done in a specific order (longest stride first, unit-stride last), and the intermediate results must be multiplied by the correct "twiddle factors" between steps.  If this is done correctly, the result will be the transpose of the pseudo-2D or pseudo-3D result, so an additional transposition is typically required.

To do a three-dimensional FFT, you simply perform a collection of 1D FFTs in each of the three directions.   If your fast 1D FFTs are set up for contiguous data, then you will need to transpose the data before each of the 1D transforms, but many packages (including MKL) are set up to handle this automatically.

0 Kudos
Masrul
Beginner
662 Views

So referring your last para of comment, I did correctly?

--Masrul

0 Kudos
McCalpinJohn
Honored Contributor III
662 Views

It looks like you have set up the DftiCreateDescriptor for a 3D transform. 

Why are you doing the "reshape" operation?   The MKL examples (e.g., https://software.intel.com/en-us/node/522290) don't do this....

0 Kudos
Masrul
Beginner
662 Views

Dr. McCalpin;

I am following  link:  https://software.intel.com/en-us/node/471390 , example number 3. Would you please look at example #3. I understand, what you are saying but above link saying otherwise.

Thanks

Masrul

0 Kudos
McCalpinJohn
Honored Contributor III
662 Views

Sorry -- I got confused by the C++ stuff.

What sort of trouble are you getting?  Is MKL complaining about the arguments or are you getting the wrong results?

0 Kudos
Masrul
Beginner
662 Views

Dr. McCalpin;

Actually i am using it  in Particle Mesh Ewald potentail energy calculation, Where i got some discrepancy at precision level or decimal level. That's why i am double checking if discrepancy from FFT part. 

--Masrul

0 Kudos
Reply