2

As the title states I'm using FFTW (version 3.2.2) with Fortran 90/95 to perform a 2D FFT of real data (actually a field of random numbers). I think the forward step is working (at least I am getting some ouput). However I wanted to check everything by doing the IFFT to see if I can re-construct the original input. Unfortunately when I call the complex to real routine, nothing happens and I obtain no error output, so I'm a bit confused. Here are some code snippets:

implicit none

include "fftw3.f"

! - im=501, jm=401, and lm=60

real*8    :: u(im,jm,lm),recov(im,jm,lm)
complex*8 :: cu(1+im/2,jm)
integer*8 :: planf,planb    
real*8    :: dv

! - Generate array of random numbers
dv=4.0
call random_number(u)
u=u*dv
recov=0.0

k=30

! - Forward step (FFT)

call dfftw_plan_dft_r2c_2d(planf,im,jm,u(:,:,k),cu,FFTW_ESTIMATE)
call dfftw_execute_dft_r2c(planf,u(:,:,k),cu)
call dfftw_destroy_plan(planf)

! - Backward step (IFFT)

call dfftw_plan_dft_c2r_2d(planb,im,jm,cu,recov(:,:,k),FFTW_ESTIMATE)
call dfftw_execute_dft_c2r(planb,cu,recov(:,:,k))
call dfftw_destroy_plan(planb)

The above forward step seems to work (r2c) but the backward step does not seem to work. I checked this by differencing the u and recov arrays - which ended up not being zero. Additionally the max and min values of the recov array were both zero, which seems to indicate that nothing was changed.

I've looked around the FFTW documentation and based my implementation on the following page http://www.fftw.org/fftw3_doc/Fortran-Examples.html#Fortran-Examples . I am wondering if the problem is related to indexing, at least that's the direction I am leaning. Anyway, if any one could offer some help, that would be wonderful!

Thanks!

JRC
  • 63
  • 2
  • 7
  • 1
    One thing to be aware of is that when you do an FFT followed by an IFFT the result will be scaled by a constant factor relative to the original input. The value of the scaling factor depends on how the FFT and IFFT are defined/implemented but it's typically N, where N is the number of points in the FFT. Of course you may have other unrelated problems here too. – Paul R Sep 22 '11 at 20:12
  • 1
    I figured out the issue - it is entirely unrelated to the FFT calls. In my code (not the above example), I had made k a 1D array but had been accessing it as a scalar index value. There shouldn't be anything wrong with the above code. Woops! Just so this is still useful for other folks, I did end up finding some code online that contained many very helpful examples of calling FFTW routines with Fortran: http://people.sc.fsu.edu/~jburkardt/f_src/fftw3/fftw3_prb.f90 – JRC Sep 22 '11 at 21:16
  • @Jen In the code above you are losing precision by mixing real*8 and complex*8 variables: complex*8 corresponds to real*4, not real*8. Use complex*16 instead. – ev-br Sep 24 '11 at 10:34
  • Jen, I actually could not even get that example code to work. Did you enable anything in the start to have it work with gfortran? He is also calling a file that doesn't exist in the current version since there is only a .f wrapper as opposed to a .f90. –  Dec 09 '11 at 07:01
  • @Thomas James, I didn't try the example with gfortran - I'm using xlf instead. Are you sure you have the correct commands/syntax for your compile steps? The .f wrapper should be fine (replace the .f90). If you would like I can send you my Makefile that I used to compile the code - just let me know. I tried pasting it here but it isn't suitable as a "comment". – JRC Dec 10 '11 at 18:41

1 Answers1

1

Not sure if this is the root of all troubles here, but the way you declare variables may be the culprit.

For most compilers (this is apparently not even a standard), Complex*8 is an old syntax for single precision: the complex variable occupies a total of 8 bytes, shared between the real and the imaginary part (4+4 bytes).

[Edit 1 following Vladimir F comment to my answer, see his link for details:] In my experience (i.e. the systems/compiler I ever used), Complex(Kind=8) corresponds to the declaration of a double precision complex number (a real and an imaginary part, both of which occupy 8 bytes).

On any system/compiler, Complex(Kind=Kind(0.d0)) should declare a double precision complex.

In short, your complex array does not have the right size. Replace occurences of Real*8 and Complex*8 by Real(kind=8) and Complex(Kind=8) (or Complex(Kind=kind(0.d0)) for a better portability), respectively.

BenBoulderite
  • 336
  • 1
  • 9