Decomposition of 2D FFT (and inverse fft) with fftw in fortran77

Question

[edit 1] Added figures to show the original data and the obtained data

[edit 2] I found my mistakes, I used fftw_measure instead of fftw_estimate in the call of dfftw_plan_many_dft

[edit 3] corrected a typo in the code (replace u with u2d in dfftw_execute_dft_r2c )

I'm trying perform the 2D fft of an array using multiple 1D fft instead of using the 2D fft function already existing in the fftw library. And subsequently, I need to perform the inverse 2D fft. The reason for doing that is that (in the future) my array will be too big to be loaded in one go to perform the 2D fft.

The 1st draft of my code looks more or less like this at the moment:

  double precision u2d(nx,ny),u2d2(nx,ny)
  double complex qhat2d(nx/2+1,ny),qhat2d2(nx/2+1,ny)
  integer N(1)
  integer howmany, idist, odist, istride, ostride
  integer inembed(2), onembed(2)
  integer rank

  ! some function to read the data into u2d

  ! perform x-fft
  N(1) = NX
  howmany = NY
  inembed(1) = NX
  inembed(2) = NY
  istride = 1
  idist = NX
  ostride = 1
  odist = (NX/2+1)
  onembed(1) = (NX/2+1)
  onembed(2) = NY
  rank = 1
  write(*,*) 'u', u2d(1,1)
  CALL dfftw_plan_many_dft_r2c(PLAN,rank,N(1),howmany,
&                              u2d,inembed,
&                              istride,idist,
&                              qhat2d,onembed,
&                              ostride,odist,FFTW_ESTIMATE) ! 
  CALL dfftw_execute_dft_r2c(PLAN,u2d,qhat2d) ! x-fft
  CALL dfftw_destroy_plan(PLAN)


  ! perform y-fft
  N(1) = NY
  howmany = (NX/2+1)
  inembed(1) = (NX/2+1)
  inembed(2) = NY
  istride = (NX/2+1)
  idist = 1
  ostride = (NX/2+1)
  odist = 1
  onembed(1) = (NX/2+1)
  onembed(2) = NY
  rank = 1
  CALL dfftw_plan_many_dft(PLAN,rank,N(1),howmany,
&                              qhat2d,inembed,
&                              istride,idist,
&                              qhat2d2,onembed,
&                              ostride,odist,FFTW_FORWARD,
&                              FFTW_MEASURE) ! 
  CALL dfftw_execute_dft(PLAN,qhat2d,qhat2d2) ! y-fft
  CALL dfftw_destroy_plan(PLAN)

  ! normally here, perform some filtering operation
  ! but at the moment, I do nothing

  ! perform inv-y-fft
  N(1) = NY
  howmany = (NX/2+1)
  inembed(1) = (NX/2+1)
  inembed(2) = NY
  istride = (NX/2+1)
  idist = 1
  ostride = (NX/2+1)
  odist = 1
  onembed(1) = (NX/2+1)
  onembed(2) = NY
  rank = 1
  CALL dfftw_plan_many_dft(PLAN,rank,N(1),howmany,
 &                             qhat2d2,inembed,
 &                              istride,idist,
 &                              qhat2d,onembed,
 &                              ostride,odist,FFTW_BACKWARD,
 &                              FFTW_MEASURE) ! 
  CALL dfftw_execute_dft(PLAN,qhat2d2,qhat2d) ! inv-y-fft
  CALL dfftw_destroy_plan(PLAN)

  ! perform inv-x-fft
  N(1) = NX ! I'm not too sure about this value here
  howmany = NY
  inembed(1) = (NX/2+1)
  inembed(2) = NY
  istride = 1
  idist = (NX/2+1)
  ostride = 1
  odist = NX
  onembed(1) = NX
  onembed(2) = NY
  rank = 1
  CALL dfftw_plan_many_dft_c2r(PLAN,rank,N(1),howmany,
&                              qhat2d,inembed,
&                              istride,idist,
&                              u2d2,onembed,
&                              ostride,odist,FFTW_ESTIMATE) ! 
  CALL dfftw_execute_dft_c2r(PLAN,qhat2d,u2d2) ! x-fft
  CALL dfftw_destroy_plan(PLAN)
  write(*,*) 'u2d2', u2d2(1,1)

  do i=1,nx
   do j=1,ny
    u2d2(i,j) = u2d2(i,j) / (nx*ny)
   enddo
  enddo
  write(*,*) 'u2d2', u2d2(1,1) ! here the values u2d2(1,1) is different from u2d(1,1)

  ! some action to write u2d2 to file
  end

I was expecting u2d and u2d2 to be the same but I obtain relatively different values. Did I make a mistake somewhere?

The original and results are shown hereunder. The shape looks similar but the values are relatively different (the minimum and maximum for example).

Original field

Field obtained after fft and i-fft

Answer 1

To clear the confusion. What happens is that the FFTW c2r and r2c routines are not guaranteed to preserve the input arrays. They can overwrite the result with garbage.

Now, you can be lucky and just use FFTW_ESTIMATE instead of FFTW_MEASURE, but it is not a good idea.FFTW_MEASURE tries many algorithms and is therefore likely to try also one which does not preserve the input. FFTW_ESTIMATE will not try to compute anything and will not overwrite the input.

The problem is that your input can be overwritten when performing the transform (executing the plan). You will only be lucky when FFT_ESTIMATE selects an algorithm which preserves the input for you. It is a lottery.

To be sure the input gets preserved you should use FFTW_INPUT_PRESERVE in addition to either FFT_ESTIMATE or FFTW_MEASURE.

You can also not use it and instead save the input somewhere. I think that is often better because FFTW_INPUT_PRESERVE can (quite likely) select a slower algorithm.

Answer 2

I found my mistakes, I used fftw_measure instead of fftw_estimate in the call of dfftw_plan_many_dft.

Correcting that gives me the appropriate original field.