FFTW fortran 90: allocatable cut in half when N>20

Question

Hi everyone and happy new year !

I'm trying to use the fftw library in a simple fortran 90 code (yes, an old fortran...). This is a very simple code computing the FFT of vector in=1,2,..., N. I'm surprise by the fact that, for N<20, it works. For N >= 20, it does not work anymore. I guess I missed something important but can't figure out what... And was wondering if you could help me... I compile my code with this command

ifort test.f90 -o test -lfftw3f

And the code is the following

program test

implicit none
include "fftw3.f"
integer, parameter :: fp =4
integer*8 :: N
double complex, allocatable, dimension (:) :: in, out, aux
integer*8 :: plan
integer*8 :: i, errflag

N=10
allocate(in(N), stat=errflag)
allocate(out(N), stat=errflag)

do i=1,N
in(i) = i
end do

call sfftw_plan_dft_1d(plan, N, in, out, -1, 0)

do i=1,N
print *, in(i)
end do
print *, "================================================"
do i=1,N
print *, out(i)
end do
call sfftw_execute_dft(plan, in, out)
call sfftw_destroy_plan(plan)

deallocate(in, out)

end program test

Surpringly (for me), the vector "in" is modified after the line

call sfftw_plan_dft_1d(plan, N, in, out, -1, 0)

Indeed, the vector is "cut in half" as soon as N>20, in the sense that:

in(i) = 0 if i < N/2
in(i) = i otherwise

However, with N =10 for example, the result seems to be good (same as the one obtained with scilab fft function).

I'm kind of lost and not totally familiar with fortran. Did I missed something important ?

Thank you so much in advance !

edit : whoups, bad copy/paste in the code...

Answer 1

Looking around for what your flags meant, I found this here http://www.fftw.org/fftw3_doc/Planner-Flags.html#Planner-Flags

Important: the planner overwrites the input array during planning unless a saved plan (see Wisdom) is available for that problem, so you should initialize your input data after creating the plan. The only exceptions to this are the FFTW_ESTIMATE and FFTW_WISDOM_ONLY flags, as mentioned below.

Maybe try

call sfftw_plan_dft_1d(plan, N, in, out, -1, 0)

do i=1,N
in(i) = i
end do

Or I may be misreading something completely out of contect, but I guess it's worth a try :)