Why would a computational scientist need to implement their own version of std::complex?

Question

Many of the better-known C++ libraries in computational science such as Eigen, Trilinos, and deal.II use the standard C++ template header library object, std::complex<>, to represent complex floating-point numbers.

In Jack Poulson's answer to a question about default constructors, he points out that he has his own implementation of std::complex in Elemental "for a number of reasons". What are those reasons? What are the advantages and disadvantages to this approach?

Answer 1

I believe this discussion has come up a number of times on the PETSc list. My main reasons are:

1. The C++ standard states that std::complex is only defined for the float, double, and long double datatypes. Thus it cannot be used for other datatypes, such as quad-precision.

2. The standard makes no guarantees about the stability of the complex arithmetic.

3. The standard makes no guarantee that the data in an std::complex is stored as the real component followed by the imaginary component. This is crucial for interfaces with external libraries, such as BLAS and LAPACK. It is true for all major implementations, but I would prefer to be able to ensure it.

4. I prefer being able to directly manipulate the real and imaginary components. std::complex makes this unnecessarily difficult.

5. I would like to eventually have a more general version which only requires the datatype to be a ring instead of requiring a field. This would include the Gaussian integers.

Answer 2

I use std::complex<> in my programs, and have to fight with compiler flags and workaround for each new compiler or compiler upgrade. I will try to recount these fights in chronological order:

1. Performance measurements showed that a a step involving only computing the square of the absolute value of a field of complex numbers took more time than a preceding FFT for gcc-4.x. Digging into the generated assembler code showed that std::norm ($|z|^2$) computed the absolute value ($|z|$) in a way avoiding overflow, and then squared the result. This problem could be fixed by the compile flag -ffast-math.
2. The intel icc compiler on linux (or linker) compiled std::arg to a non-opt under certain configurations (link compatibility with a specific gcc-version). The problem resurfaced too often, so std::arg had to be replaced by atan2(imag(),real()). But it was all too easy to forget this when writing new code.
3. The type std::complex uses different call conventions (=ABI) than the build-in C99 complex type, and the built-in Fortran complex type for newer gcc versions.
4. The -ffast-math compile flag interacts with the handling of floating point exceptions in unexpected ways. What happens is that the compiler pulls divisions out of loops, thereby causing division by zero exceptions at runtime. These exceptions would have never happened inside the loop, because the corresponding division didn't take place due to the surrounding logic. That one was really bad, because it was a library that was compiled separately from the program which used the floating point exception handing (using different compile flags) and run into these issues (the corresponding teams were sitting at opposite parts of the world, so this issue really caused bad trouble). This was solved by doing the optimization used by the compiler by hand with more care.
5. The library became part of the program and no longer used the -ffast-math compile flag. After an upgrade to a newer gcc version, performance dropped by a huge factor. I haven't investigated this issue in detail yet, but I fear it is related to C99 Annex G. I have to admit that I'm completely confused by this strange definition of multiplication for complex numbers, and there even seems to exist different versions of this with claims that the other versions are misguided. I hope that the -fcx-limited-range compile flag will solve the issue, because there seems to be another problem related to -ffast-math for this newer gcc version.
6. The -ffast-math compile flag makes the behavior of NaN completely unpredictable for newer versions of gcc (even isnan is affected). The only workaround seems to be to avoid any occurrence of NaN in the program, which defeats the purpose for the existence of NaN.

Now you may ask whether I plan to abandon the built-in complex types and std::complex for these reasons. I will stay with the built-in types, as long as I stay with C++. In case C++ should manage to become completely unusable for scientific computing, I would rather considering to switch to a language that takes more care of the issues relevant to scientific computing.