Dear Ying,
On Thu, Jan 14, 2010 at 1:00 PM, Ying Zou
Dear all,
One colleague of mine brought up this to my attention. He tried a simple Fourier transformation C code on different machines (PowerPC based MAC G4, Intel based MAC G5, and PCs). The computations yield different results, sometimes on the third significant place. In his Fourier transformation implementation, a standard complex function is called. On MAC boxes, he used GNU C compiler; on PC, MS visual C compiler was used. I have done a little web survey on this issue, it seems that this was reported by others, too; there have been suggestions that it is much more reliable to handle real or imaginary part separately.
By "more reliable" do you mean in general or for these particular configurations (compilers, processors)? You didn't mention which libraries were used for the FTs or what precision level was expected. It does seem like the work is restricted to C implementations. You also don't give any references for the "suggestions of more reliability", so it's a bit hard to discuss those. It is possible that one of the problems is in expecting "simple" and "fourier transformation" to go together (they usually don't) or that the FT algorithm used was not robust.
Since in the Fourier transformation for EXAFS function, complex functions have to be handled too. One thing I would like to know from our Ifeffit Gurus: Are the real and imaginary parts of EXAFS functions been handled separately?
The Fourier transforms used in Ifeffit are on complex numbers. Ifeffit is written in Fortran and uses Fortran's intrinsic support for complex numbers, adhering closely to the F77 standard so that it can be compiled by g77 (in one important aspect for this discussion the strict F77 is not followed -- more later). Ifeffit uses a subset of the Fortran implementation in FFTPACK (which is included with the Ifeffit source and compiled and statically linked into Ifeffit) for discrete Fourier transforms. Arrays to be Fourier Transformed in Ifeffit are put onto a grid of 2048 points, and zero-padded as needed. These arrays are put into arrays of complex*16 -- double-precision complex (this is a non-standard extension of F77 but implemented for almost all existing compilers) which are then transformed by FFTPACKs cfftf or cfftb routines depending on whether a forward or reverse transform is desired. See http://www.netlib.org/fftpack/doc for a little more detail on FFTPACK. Although Ifeffit uses double precision and double precision complex, it does not ensure that the data or constants used as normalization factors are kept at full double precision (XAFS data is rarely measured or stored beyond single precision). I would not be surprised at all if results varied from machine to machine at single precision levels (1e.-7 in relative precision). If it was only to 1.e-4, I would not be horrified. High precision and accuracy in floating point computation is truly hard, and machine-to-machine variations are common. The key need for Ifeffit is to be able to compare the Fourier transform of two spectra computed under the same conditions. Of course, there may be bugs or ways to improve problems in precision. If you have suggestions, I'd be happy to hear them. Hope that helps. --Matt