Bug in complex numbers

Comments (13)

1. 

   Frank Löffler

   • removed comment

   From what I read, I assume when you say

   I hear that this makes Cactus complex numbers unusable on 32-bit Intel systems (untested). you mean the current implementation, right? The problem seems to be that the C++ complex<double> is not (necessarily) ABI compatible to C's "complex double" - which is the same as Fortan's?

   • 2014-02-06T14:07:56+00:00
2. 

   Frank Löffler

   • removed comment

   And just realizing this - it would also mean that the current complex number implementation is fine, if you don't use it from C++. Which of course should be corrected - just trying to understand the implications.

   • 2014-02-06T14:10:29+00:00
3. 

   Erik Schnetter reporter

   • removed comment

   I assume (but I don't know) that Fortran and C use the same ABI.

   Yes, only C++ is broken, since Cactus implements the operations in C.

   • 2014-02-06T14:16:01+00:00
4. 

   Roland Haas

   • removed comment

   Reading the bug report this seems to be only an issue for return values. The bug report states that even std::complex<double> is fine in memory ie std::complex<double> must look the same as _Complex double in memory. The only functions we have that returns a CCTK_COMPLEX are the CCTK_CmplxAdd etc functions. Since we now require C99 all of our supported languages (C, C++, F77, F90) support a complex type natively to do algebra with. So it seems feasible that we can make things standards conforming by removing the CCTK_CmplxAdd etc. functions which are no longer required since the users can just use the language provided functionality instead. For backwards compatibility we can provide macro versions of the routines in cctk_Math.h or similar (or inline functions returning CCTK_COMPLEX for C/C++). I actually do run the Cactus test suites on a 32bit (virtual) machine. The failing tests (current trunk) are (using gcc/gfortran):

       waveinterp-2p (from CarpetInterp)
       test_o9 (from CarpetProlongateTest)
       Schwarzschild_EF (from Exact)
       test_tov_carpet_refined_nosync (from GRHydro)
       tov_slowsector (from GRHydro)
       test_axibrill_nostagger (from IDAxiBrillBH)
       test_simpson (from Multipole)
       regression_test (from SphericalHarmonicRecon)
   

   SphericalHarmonicRecon in PITTNull uses complex numbers. PITTNull itself is in Fortran so probably fine, MoL does the update and is in C so should be fine as well. There seems to be no C++ code using complex numbers in the tests. Slab uses complex numbers but the test does not use Slab. Both SphericalHarmonicReconGen and SphericalSlice do use them and are in C++, neither one is in the ET yet.

   There are two reasons I can see for the SphericalHarmonicRecon test failing: 1. complex return values are handled differently in C and Fortran after all (different from what is stated in the link above) 2. 32bit machines by default compile for the i387 coprocessor which always uses 80bit long double precision in its registers (well almost always) while 64bit machines use sse2 which uses only 64bit double precision it its registers. One can test this by recompiling the code (assuming gcc) using -mfpmath=sse

   I am testing patch for the flesh that replaces the complex algebra routines by inline function calls which should remove the problem of possibly differing return value representations and will report back on the results using the patch alone and using -mfpmath=sse.

   • 2014-02-06T17:56:15+00:00
5. 

   Erik Schnetter reporter

   • removed comment

   I checked: On a 32-bit system, C and Fortran use the same calling convention, and C++ differs. Thus it seems we can't use C++ complex numbers.

   We currently typedef CCTK_COMPLEX differently in C and C++. I believe we need to change this, using "_Complex double" in both languages. Otherwise, e.g. a user-defined C routine returning CCTK_COMPLEX would fail when called from C++ and vice versa.

   The best solution I can come up with is to introduce routines to the flesh that convert from/to CCTK_COMPLEX and complex<CCTK_REAL>. User code would then have to call these explicitly when accessing grid functions.

   We can also introduce a new type CCTK_COMPLEX_CXX (namely complex<CCTK_REAL>), and use it to declare complex grid functions in C++. We can define this type even for C (as a struct, which has the correct memory layout). However, this would probably lead to trouble with header files included from both C and C++. I don't think this is a good idea.

   As mentioned above, the standard operators should work in C++ even for C complex types (_Complex double). We could provide our own overloads of std::abs for these, so that C++ code "would not notice". Maybe these overloads should be in the cctk namespace instead.

   • 2014-02-06T18:40:23+00:00
6. 

   Roland Haas

   • changed status to open
   • removed comment

   On comment:5: introducing a struct type for CCTK_COMPLEX in C seem undesirable. We'd loose the ability to to arithmetic on them even when using C99. We'd also have no guarantee anymore that it will look the same as Fortran's complex numbers (though I'd be surprised if they would), nor does C++ seem to require that the compiler actually implements complex numbers as a user visible class. As far as I remeber C++ allows the compiler to recognize "#include <complex>" and then map std::complex<double> to some internal representation about which is does not have to tell us anything.

   I checked the test results and the test failures in the ET when using 32bit machines are indeed due to the 80bit math precision. Compiling with -mfpmath=sse -march=prescott (and -O0 though I doubt that matters since it alone did not help) I get the test (in PITTNull) to pass.

   Attached please find a patch for the flesh that removes the complex algebra helper routines and replaces them by inlined copies, one for C++ and one for C. This leaves the case where a user defines C routine that returns CCTK_COMPLEX but calls it from C++ where we I don't know which calling convention is used for a routine declared extern "C".

   • 2014-02-06T19:10:22+00:00
7. 

   Erik Schnetter reporter

   • removed comment

   I attach an updated version of the patch that does not include <complex.h> from header files.

   • 2014-02-07T14:17:12+00:00
8. 

   Roland Haas

   • removed comment

   Are the changes to Complex.c (compared to the original trunk version) required? It seems as if the new code does the same as the code it replaces only with fewer comments. The remove all the grdocs text which (in the flesh) are in our coding guide. If sufficient, I'd only change cctk_Complex.h to provide inline functions for C++ and the comment about not pushing "I" into global address space in front of the C versions.

   • 2014-02-07T15:28:09+00:00
9. 

   Erik Schnetter reporter

   • removed comment

   You are right; the changes to Complex.c should not be necessary.

   • 2014-02-07T15:41:27+00:00
10. 

    Erik Schnetter reporter

    • removed comment

    Do you agree with the complex.2.patch if Complex.c remains unchanged?

    • 2014-02-08T13:40:30+00:00
11. 

    Roland Haas

    • changed status to open
    • removed comment

    Yes I agree with complex.2.patch if Complex.c remains unchanged. I attached such a version as complex.3.patch. Passes the tests on my 32bit virtual machine.

    Please apply.

    • 2014-02-08T15:13:38+00:00
12. 

    Erik Schnetter reporter

    • changed status to resolved
    • removed comment

    Applied.

    • 2014-02-11T15:27:25+00:00
13. 

    Roland Haas

    • edited description
    • changed status to closed

    • 2019-02-21T20:11:32+00:00
14. Log in to comment