2D FFT on 4D data real to complex and backward

[the_t_1]

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

[Chao_Y_Intel]

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

[Dmitry_B_Intel]

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

[the_t_1]

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

[the_t_1]

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.

[Dmitry_B_Intel]

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.