- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
Please can I ask for help with the 2D transform on 4D data? I have an array theta(x,y,dz,z) and need to do to the FFT over x and y, I have written for this an test code:
real :: Xin(x,y), Xout(x*y)
complex :: Yin(x,y), Yout(x*y)
dxT = x/2 + 1
L(1) = x
L(2) = y
strides_out(1) = 0
strides_out(2) = 1
strides_out(3) = dxT
do i=1,dz
do j = 1,z
Xin = theta(:,:,i,j)
StatusExp = DftiCreateDescriptor( FFT_HANDLE, DFTI_SINGLE,DFTI_REAL, 2, L )
StatusExp = DftiSetValue(FFT_HANDLE,DFTI_CONJUGATE_EVEN_STORAGE, DFTI_COMPLEX_COMPLEX)
StatusExp = DftiSetValue( FFT_HANDLE, DFTI_PLACEMENT, DFTI_NOT_INPLACE )
StatusExp = DftiSetValue(FFT_HANDLE, DFTI_OUTPUT_STRIDES, strides_out)
StatusExp = DftiCommitDescriptor(FFT_HANDLE)
StatusExp = DftiComputeForward(FFT_HANDLE, Xin(:,1), Yout)
do k=1,dxT
do l=1,y
thetaFFT(k,l) = Yout(k + (l-1)*dxT)
if(k>1) thetaFFT(x+2-k,l) = conjg(thetaFFT(k,l))
enddo
enddo
StatusExp = DftiFreeDescriptor(FFT_HANDLE)
enddo
enddo
!***************************************************************************
! Backward
L(1) = x
L(2) = y
strides_out(1) = 0
strides_out(2) = 1
strides_out(3) = x
do i=1,dz
do j = 1,z
Yin = thetaFFT(:,:,i,j)
StatusExp = DftiCreateDescriptor( FFT_HANDLE, DFTI_SINGLE,DFTI_REAL, 2, L )
StatusExp = DftiSetValue(FFT_HANDLE,DFTI_REAL_STORAGE, DFTI_REAL_REAL)
StatusExp = DftiSetValue( FFT_HANDLE, DFTI_PLACEMENT, DFTI_NOT_INPLACE )
StatusExp = DftiSetValue(FFT_HANDLE, DFTI_OUTPUT_STRIDES, strides_out)
StatusExp = DftiCommitDescriptor(FFT_HANDLE)
StatusExp = DftiComputeBackward(FFT_HANDLE, Yin(:,1), Xout)
do k=1,x
do l=1,y
theta(k,l) = Xout(ix + (l-1)*dxT)/(x*y)
enddo
enddo
StatusExp = DftiFreeDescriptor(FFT_HANDLE)
enddo
enddo
Why I don't get my testing array theta back when running this routine?
I send theta = 1., procedure with it forward and backward FFt and get back a theta containing random numbers about 2. What's wrong?
Many thanks for any idea
Link Copied
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
Hello,
For backward transform, you are setting the DFTI_REAL_REAL, DFTI_REAL_REAL is only used store the real data:
DftiSetValue(FFT_HANDLE,DFTI_REAL_STORAGE, DFTI_REAL_REAL)
You can use DFTI_CONJUGATE_EVEN_STORAGE, storage, Yin set the same data as the Yout:
You can check this article for some example:
http://software.intel.com/sites/products/documentation/hpc/mkl/mklman/GUID-0E5501C5-30C5-413F-BD13-FC45CBC2884A.ht
Thanks,
Chao
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
Hi the t,
Conjugate-even symmetry of the result of the first transform should be employed by this:
thetaFFT(x+2-k,y+2-l) = conjg(thetaFFT(k,l))
Thanks
Dima
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
Many thanks, Chao and Dima,
I have implemented your advices but it still didn't work. So I have written a very simple code and have found, that when I compute tha Forward and Backward in one loop, it works. But not, when I close the loop and open new! But I need to compute the whole theta(x,y,dz,z) and do some algebra with it and then finally convert it back to direct space. So I need to close the loops..
Here is the program which works properly:
! first, theta is filled with random numbers, symbolically:
! theta = rand()
L(1) = dx
L(2) = dy
strides_out(1) = 0
strides_out(2) = 1
strides_out(3) = dx
do i=1,idz
do j=1,iz
xy_data = cmplx(theta(:,:,i,j),0)
StatusExp = DftiCreateDescriptor(hand, DFTI_DOUBLE, DFTI_COMPLEX, 2, L)
StatusExp = DftiCommitDescriptor(hand)
StatusExp = DftiSetValue( hand, DFTI_PLACEMENT, DFTI_NOT_INPLACE )
StatusExp = DftiSetValue( hand, DFTI_OUTPUT_STRIDES, strides_out)
StatusExp = DftiComputeForward(hand, xy_data(:,1), lambda_data(:,1))
thetaFFT(:,:,i,j) = lambda_data
StatusExp = DftiFreeDescriptor(hand)
! enddo
!enddo
!do i=1,idz
! do j=1,iz
lambda_data = thetaFFT(:,:,i,j)
StatusExp = DftiCreateDescriptor(hand, DFTI_DOUBLE, DFTI_COMPLEX, 2, L)
StatusExp = DftiCommitDescriptor(hand)
StatusExp = DftiSetValue( hand, DFTI_PLACEMENT, DFTI_NOT_INPLACE )
StatusExp = DftiSetValue( hand, DFTI_OUTPUT_STRIDES, strides_out)
StatusExp = DftiComputeBackward(hand, lambda_data(:,1), xy_data(:,1))
theta_new(:,:,i,j) = real(xy_data)
StatusExp = DftiFreeDescriptor(hand)
enddo
enddo
But if I uncomment the part in the middle,
do i=1,idz
do j=1,iz
xy_data = cmplx(theta(:,:,i,j),0)
StatusExp = DftiCreateDescriptor(hand, DFTI_DOUBLE, DFTI_COMPLEX, 2, L)
StatusExp = DftiCommitDescriptor(hand)
StatusExp = DftiSetValue( hand, DFTI_PLACEMENT, DFTI_NOT_INPLACE )
StatusExp = DftiSetValue( hand, DFTI_OUTPUT_STRIDES, strides_out)
StatusExp = DftiComputeForward(hand, xy_data(:,1), lambda_data(:,1))
thetaFFT(:,:,i,j) = lambda_data
StatusExp = DftiFreeDescriptor(hand)
enddo
enddo
do i=1,idz
do j=1,iz
lambda_data = thetaFFT(:,:,i,j)
StatusExp = DftiCreateDescriptor(hand, DFTI_DOUBLE, DFTI_COMPLEX, 2, L)
StatusExp = DftiCommitDescriptor(hand)
StatusExp = DftiSetValue( hand, DFTI_PLACEMENT, DFTI_NOT_INPLACE )
StatusExp = DftiSetValue( hand, DFTI_OUTPUT_STRIDES, strides_out)
StatusExp = DftiComputeBackward(hand, lambda_data(:,1), xy_data(:,1))
theta_new(:,:,i,j) = real(xy_data)
StatusExp = DftiFreeDescriptor(hand)
enddo
enddo
This code doesn't work (the final theta_new is not equal to input theta). I fill the theta field with random numbers. Hovewer, if I fill it with all 1, it works...
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
I have tried now to do the same with the example codes in MKL directory (basic_dp_real_dft_3d.f90) and have realized that changing the
call init... to
i = time()
call seed(i)
do i=1, N1
do j=1, N2
do k=1, N3
x_real(i,j,k) = rand()
enddo
enddo
enddo
does the same like my own code - the results are not the same. Maybe is this a property of function rand()? But I don't understand, how this may be explained.
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
The 'program which works properly' is not supposed to do so, because it changes the descriptor after it has been committed.
The sequence of steps when using DFTI is this: create a descriptor, optionally set configuration in the descriptor, commit the descriptor and check the status returned by the commit function. Then the descriptor can be used to compute the FFT arbitrary number of times. Finally free the descriptor. Please refer to documentation.

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