Hi,
I spend quite a long time finding why in fortran my FFT subroutine FOUR1 is producing wrong results, and I've seen many people had problems with that, so I would like to share my experience.


In my quite older Numerical Recipes in Fortran (1989) in FOUR1 there is a declaration of field DATA as : DIMENSION DATA(2*NN)
This was producing on my computer obviously wrong results. Good example is (as I found here in thread 'four1 output'):
for NN=4, the real sequence 0,1,2,3
so input
x(1)=0
x(2)=0
x(3)=1
x(4)=0
x(5)=2
x(6)=0
x(7)=3
x(8)=0
because complex parts are zero. This produced wrong output

1.000000
-1.000000
1.000000
-1.000000
2.000000
0.000000
3.000000
0.000000

But when declaration of the field was changed to double precision
REAL(KIND=8) :: DATA(2*NN)
or
DOUBLE PRECISION DATA(2*NN)
the results become correct (as verified with matlab and fftpack), with output

6.000000
0.000000
-2.000000
-2.000000
-2.000000
0.000000
-2.000000
2.000000

I'm using Intel Fortran compiler...
Same bad results as before are obtained if KIND=4 in declaration of DATA is used, so my suspicion that just my compiler is missinterpreting a little old fashioned statemenent 'dimension' is not true.

Explanation:

Then I found the reason, what many people, including me, overlooked:
Everyting depends how you in the main program declare your field what you are using for data input. When your input field is let's say X, so in the main program you are calling let's say CALL FOUR1(X,N,1), if field X has the same declaration of precision as DATA in FOUR1, you will obtain correct results. But as was in my case, I declared field X (as always) for kind=8 and FOUR1 with statement 'DIMENSION' is assuming for DATA kind=4. It is kind of strange why this declaration for data was choosen, because otherwise all FOUR1 is written everywhere with double precision...

CONCLUSION:

For data, use the same kind=4 or kind=8 both in main program and subroutine FOUR1 and your results will be correct.

Maybe this is was for many people obvious, but I almost threw my computer out of window... :)