Re: What is acceptable for -ffast-math? A numerical viewpoint

[Wolfgang Bangerth <wolfgang dot bangerth at iwr dot uni-heidelberg dot de>]

Due to repeated demand: here is the opinion of a numerical analyst. To
further define my position: there is some kind of schism in the numerics
society, the "traditionalists" doing Fortran with highly tuned codes and
relatively simple data structures, and the "modernists" using extremely
large C++ programs on complex data structures. I belong to the latter
group and can only speak for them. I know that the former group has
different views on these topics.

My opinion is that -fast-math could well include slight deviations from
IEEE modes for denormals, rounding modes, associative redistribution, etc.
The reasons being:
- denormals: if you run into the risk of getting NaNs one way but not
  the other, you most certainly already have gone beyond the limits of
  your program, i.e. your linear solver diverged, you have chosen an
  inherently unstable formulation for your problem, etc. In practice, if
  you don't get you NaN on this operation, you will get it some operations
  later, and if not then the result is spoilt anyway because of the
  instability. I cannot imagine someone writing a code which runs on
  a well-defined side of the boundary between normal and denormalized
  values and still returns reasonable results. And even if such a program
  existed, then that person must have known exactly what he does and will
  refrain from -fast-math and use -mieee instead.
- rounding modes: if it makes a difference, then again you are already
  in an instable regime. Then why bother: I know that I only solve an
  approximation of the true laws of nature (e.g. after discretizing a
  partial differential equation), why solve that to 16 digits of accuracy?
- association: Fortran codes differ significantly from C++ codes in
  some ways, not the least are that C++ codes are often split into many
  very small parts. In these cases, I have no control over how the
  operands associate, if they come from different functions. For example,
  if I have (a very much simplified situation)
    double x, y;
    double f () { return x; };
    double g () { return y; };
    double h () { return x; };
    ...
    double z = f() * g() / h();
  then I'd be happy if I get z==y, although that might not be right
  from an IEEE viewpoint, or if f()+g()+h() would be turned into
  2*x+y. If you think that this example is unrealistic,
  think about using the C++ valarray<> class with its small functions
  combined with a case where some of the function's arguments are constant
  and known so that most of the operations in the body of the functions
  can be optimized away. If the compiler can help optimize such cases
  across the boundaries of single functions, then that would be great, but
  it will often require violating associative laws.

  Note that this is a case where C++ differs significantly from Fortran:
  there, you usually have large functions and optimizations take place
  within these.

So, concluding: if you have programs that run for several days, you'd be
happy if you could cut that time by some hours using optimizations that
- still do what common sense would dictate (IEEE is not pure common sense,
  but rather a definition), i.e. swapping operands would be allowed, but
  replacing NaNs by zero would not
- can be switched off in case I know I have to treat boundary cases
  exactly.
I myself would certainly appreciate a compiler turning a*c+b*c into
(a+b)*c, I'd definitely not be shocked about that.

Regards
  Wolfgang


-------------------------------------------------------------------------
Wolfgang Bangerth          email: wolfgang.bangerth@iwr.uni-heidelberg.de
                             www: http://gaia.iwr.uni-heidelberg.de/~wolf