[Ifeffit] Compiler cannot handle complex function properly

[Matt Newville]

Dear Ying,


On Thu, Jan 14, 2010 at 1:00 PM, Ying Zou <zou at uwm.edu> wrote:
> 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