Floating-point optimizations - guideline

Question

The majority of scientific computing problems that we need solve by implementing a particular algorithm in C/C++ demands accuracy that are much lower than double precision. For example, 1e-6, 1e-7 accuracy covers 99% of the cases for ODE solvers or numerical integration. Even in the rare cases when we do need higher accuracy, usually the numerical method itself fail before we can dream of reaching an accuracy that is near double precision. Example: we can't expect 1e-16 accuracy from a simple Runge–Kutta method even when solving a standard nostiff ordinary differential equation because of roundoff errors. In this case, the double precision requirement is analogous of as asking to have a better approximation of the wrong answer.

Then, aggressive float point optimizations seems to be a win-win situation in most cases because it makes your code faster (a lot faster!) and it does not affect the target accuracy of your particular problem. That said, it seems remarkable difficult to make sure that a particular implementation/code is stable against fp optimizations. Classical (and somewhat disturbing) example: GSL, the GNU scientific library, is not only the standard numerical library in the market but also it is a very well written library (I can't imagine myself doing a better job). However, GSL is not stable against fp optimizations. In fact, if you compile GSL with intel compiler, for example, then its internal tests will fail unless you turn on -fp-model strict flag which turn off fp optimizations.

Thus, my question is: are there general guidelines for writing code that is stable against aggressive floating point optimizations. Are these guidelines language (compiler) specific. If so, what are the C/C++ (gcc/icc) best practices?

Note 1: This question is not asking what are the fp optimizations flags in gcc/icc.

Note 2: This question is not asking about general guidelines for C/C++ optimization (like don't use virtual functions for small functions that are called a lot).

Note 3: This question is not asking the list of most standard fp optimizations (like x/x -> 1).

Note 4: I strongly believe this is NOT a subjective/off-topic question similar to the classical "The Coolest Server Names". If you disagree (because I am not providing a concrete example/code/problem), please flag it as community wiki. I am much more interested in the answer than gaining a few status points (not they are not important - you get the point!).

Rossi

Errors accumulate - even if all calculations at done at double precision the final result will not be accurate to the last bit - far from it. If you use float everywhere you'd need to do a proper error analysis to work out how many bits of your answer are trustworthy (if any). You should do the same thing for double, of course. – Thaxter

General numerical stability is often achieved via carefully chosen and mostly fragile intermediate steps, specifically designed to overcome the non-associative character of the finite-precision floating-point arithmetic. Aggressive optimisations tend to break that e.g. by changing the order of execution. Iterative adjustment could help, but then it takes more time to get the answer. Pro tip: you might want to post your question on the Computational Science site. – Selfregard

Single precision with -ffast-math works just fine for me for all purposes I have come across so far. If you choose a good algorithm and write the software well it should work. – Graniela

Answer 1

Compiler makers justify the -ffast-math kind of optimizations with the assertion that these optimizations' influence over numerically stable algorithms is minimal.

Therefore, if you want to write code that is robust against these optimizations, a sufficient condition is to write only numerically stable code.

Now your question may be, “How do I write numerically stable code?”. This is where your question may be a bit broad: there are entire books dedicated to the subject. The Wikipedia page I already linked to has a good example, and here is another good one. I could not recommend a book in particular, this is not my area of expertise.

Note 1: Numerical stability's desirability goes beyond compiler optimization. If you have choice, write numerically stable code even if you do not plan to use -ffast-math-style optimizations. Numerically unstable code may provide wrong results even when compiled with strict IEEE 754 floating-point semantics.

Note 2: you cannot expect external libraries to work when compiled with -ffast-math-style flags. These libraries, written by floating-point experts, may need to play subtle tricks with the properties of IEEE 754 computations. This kind of trick may be broken by -ffast-math optimizations, but they improve performance more than you could expect the compiler to even if you let it. For floating-point computations, expert with domain knowledge beats compiler every time. On example amongst many is the triple-double implementation found in CRlibm. This code breaks if it is not compiled with strict IEEE 754 semantics. Another, more elementary algorithm that compiler optimizations break is Kahan summation: when compiled with unsafe optimizations, c = (t - sum) - y is optimized to c = 0. This, of course, defeats the purpose of the algorithm completely.

Divebomb

It is not impossible for a library to work with "-ffast-math". One of the milestone targets of the GSL developers is to make GSL stable under -ffast-math (according to a report on gsl mail list). They say it would be hard (but not impossible). In addition, GSL only implements very reliable and stable algorithms. So it is a code issue that is not impossible to solve => there must be general techniques in C/C++ to achieve that. – Rossi

@ViniciusMiranda I said “you cannot expect external libraries”. It is a general statement about all floating-point libraries. If you prefer: Not Forall floating-point library F, you can expect F to work with -ffast-math. If the authors of GSL aim for compatibility with these options, good for them. Not everyone does and in some cases, it does not even make sense, for the kind of reason I alluded to. – Divebomb