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)
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.
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....
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.
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.