ochubar / SRW

Synchrotron Radiation Workshop
Other
86 stars 70 forks source link

fftw srw bug? #9

Closed buzmakov closed 7 years ago

buzmakov commented 8 years ago

This pull request should fix fftw bug (https://github.com/buzmakov/SRW/issues/4) discovered by @SergeyYakubov Could you please merge this changes in SRW code?

sergeyyakubov: in function SetRepresFT (srradstr.cpp) we have:

https://github.com/ochubar/SRW/blob/master/cpp/src/core/srradstr.cpp#L3640

.... if(pBaseRadX != 0) { FFT1DInfo.pInData = pBaseRadX; FFT1DInfo.pOutData = pBaseRadX; if(result = FFT1D.Make1DFFT(FFT1DInfo)) return result; } ....

inside Make1DFFT fftw is called:

https://github.com/ochubar/SRW/blob/master/cpp/src/ext/genmath/gmfft.cpp#L374

fftw(Plan1DFFT, FFT1DInfo.HowMany, DataToFFT, 1, Nx, OutDataFFT, 1, Nx);

In this case DataToFFT==OutDataFFT and FFTW_IN_PLACE flag is set. But having OutDataFFT not NULL when using in-place fftw means that OutDataFFT should be an auxilary array used:

In-place transforms: If the plan specifies an in-place transform, ostride and odist are always ignored. If out is NULL, out is ignored, too. Otherwise, out is interpreted as a pointer to an array of n complex numbers, that FFTW will use as temporary space to perform the in-place computation. out is used as scratch space and its contents destroyed. In this case, out must be an ordinary array whose elements are contiguous in memory (no striding). 

Therefore OutDataFFT (which points to DataToFFT) will be used as scratch space. This looks wrong to me. The correct call in this case should be

fftw(Plan1DFFT, FFT1DInfo.HowMany, DataToFFT, 1, Nx, 0, 0, 0);

I get slightly (which is a matter of luck) different results when using original version vs corrected.