16

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?

Tyberius
  • 1,023
  • 1
  • 10
  • 26
Aron Ahmadia
  • 6,951
  • 4
  • 34
  • 54

3 Answers3

19

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.

Jack Poulson
  • 7,599
  • 32
  • 40
  • 7
    Point 3 has been addressed in C++11. 26.4.4 states that if z is an lvalue expression of type cv std::complex<T> then reinterpret_cast<cv T(&)[2]>(z) and reinterpret_cast<cv T(&)[2]>(z)[0] shall designate the real part of z, and reinterpret_cast<cv T(&)[2]>(z)[1] shall designate the imaginary part of z. Arrays of complex numbers are also addressed. – James Custer Apr 21 '12 at 16:07
  • 3
    @JamesCuster: I'm all for eventually switching to C++11, but scientific codes that want to remain portable to semi-exotic architectures will probably need to wait at least two to three more years to do so. Also, C++11 unfortunately only addresses part of the problem. – Jack Poulson Apr 21 '12 at 21:01
  • I understand, I was just throwing it out there in case someone looks at this question in the future. – James Custer Apr 21 '12 at 21:05
  • 2
    Well, I think it's a cop-out to say you'd have to wait till compilers support C++11. The explicit requirement was put into the new standard because all existing implementations already support it. I can't think of a case where it would be unsafe to already assume this particular layout in existing compilers/libraries as it would simply not have made any sense to implement std::complex in any other way. – Wolfgang Bangerth Apr 22 '12 at 02:43
  • 1
    @WolfgangBangerth: It was more of a general comment on switching to C++11. Either way, C++11 does not fix most of the problems with std::complex. – Jack Poulson Apr 22 '12 at 05:17
  • Point 2: I just looked up the implementation of the multiplication in gcc: https://github.com/gcc-mirror/gcc/blob/9d8e27fac3c6426fbfaf30d51cbf2761c5491a6a/libstdc%2B%2B-v3/include/std/complex#L301 ; it's the "grade school" implementation. IMO I would use Kahan's difference of products algorithm to get 1.5 ulps, but that's getting into a design decision which arguably should be solved differently in different contexts. – user14717 Jan 13 '22 at 16:53
8

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.

Thomas Klimpel
  • 2,141
  • 15
  • 35
  • Looks like my fears related to C99 Annex G have come true, and -fcx-limited-range is now sort of required for decent computation speed when multiplying complex numbers. At least that is what I get from the following recent war story: https://medium.com/@smcallis_71148/complex-arithmetic-is-complicated-873ec0c69fc5 – Thomas Klimpel Aug 02 '17 at 11:48
1

Another possible reason might be a historical one: your own version pre-dates the current accepted standard.

While this isn't about complex numbers, this discussion of a bug report is on the use of a container-type. The software in question, i.e., OpenFOAM, has been developed from a time long before the standard template library (STL) of C++ has been a thing. Hence, this software contains a lot of data types, which would nowadays simply be taken from the STL.

While the sunken cost fallacy is usually cited as being a fallacy, there is nevertheless a sunken cost, which causes many developers to maintain their own implementations of "standard" types, because a refactoring of their code would take a lot of resources.

Dohn Joe
  • 559
  • 2
  • 9