Incorrect Idamax with NaN values
btracey opened this issue · 9 comments
In the reference implementation, the arrays [5, +inf, NaN, 8, 9] and [5, NaN, +inf, 8, 9] have an Idamax of 2 and 3 respectively. OpenBLAS returns 2 and 2. While this doesn't strictly have to be a bug (since the behavior with NaN isn't specified), it seems like x[idamax(x)] shouldn't depend on slice order.
NaN is not comparable. You should get floating point exception on NaN comparison and undefined result. correct answer is a guess from list of all NaNs in array and idamax of the rest. If you read idamax documentation...
Idamax documentation (http://www.netlib.org/lapack/explore-html/dd/de0/idamax_8f_source.html): "IDAMAX finds the index of the first element having maximum absolute value." It does not specify behavior when NaN is present.
With NaN not being comparable it can be bigger than infinity. Why are you insisting on consistency on undefined behaviour?
Especially the order of comparison just cannot happen with parallel and vectorized BLAS.
inf.gt.NaN gives different results on different compilers and different flags (like gfortran fast-math vs i387 vs SSE2 vs intel vs pgi) - do you expect most popular case is correct or idamax of version numbers should be winning case?
Solution for your expectation is to wrap the idamax function and replace NaNs with zeroes.
As it stands today idamax has no chance to return error conditions like seeing nans
in real life you need to rewind your computation back to origins and NOT create NaNs first hand.
@brada4 , thank you for the update. So far, OpenBLAS is very hard to deal with NaN.
Implementation compatibility may make little sense when it comes to undefined behaviour, so how about maximum execution speed. Anybody know the most efficient way to identify NaN values in assembly (for x86 I am led to believe it is ucomisd) ? Even if the OpenBLAS *amax implementations elected to return the highest possible index number immediately upon encountering the first NaN element, it would solve #671, #723 without being "wrong". The alternative would appear to be to clearly document that OpenBLAS intentionally relies on numerically valid input.
In case of NaN-s in (netlib) *max and very lucky compiler and fpu combination it can happen NaN >= Inf, and in the next step 1e-99 >= NaN, and it is exactly what undefined behaviour means.
Sure it can be optimized to common sense that first infinity or Nan is maximum and dont search further, but then again bug report implies NaNs should be replaced with fake zeroes.
NaN check before LAPACK calls is only (optionally) present in LAPACKE, and not before that. Which means that silent garbage in garbage out is very well within specification.
GIGO is one thing, crashing is another. My suggestion is either strive for the fastest non-crashing implementation no matter if its result in the NaN case happens to coincide with netlib-du-jour, or document that OpenBLAS trades NaN safety for speed.
This one is not crashing, it is just confused by NaN in the way different from reporter's fortran compiler.
"BLAS predates IEEE 754 and does not have any provisions for concept of NaN in most cases"