DenseVector::infinity_norm() does not handle Infinity properly
I hate to bring this up again but I'm afraid I'll have to. I'll quickly lay out the history of this issue:
If you compute something like DenseVector::two_norm()
for a vector that contains a value that is not finite (either Infinity
or NaN
), that will show up in the final result. Because all the operations involved in the computation of DenseVector::two_norm()
, namely operator+=
, operator*
and std::sqrt
will pass these values through.
If instead you compute something like DenseVector::infinity_norm()
at some point you'll find yourself calling a function that computes the maximum of two numbers. That could be std::max
or std::fmax
. In some cases it makes sense to treat entries in a vector or matrix as missing, so that that they should be dropped from a computation. A language like R
has NA
for that; it's different from NaN
. With doubles in C++
, we don't have NA
, just NaN
; which is maybe why it makes sense that std::fmax
ignores NaN
arguments: It would be the logical thing to do if NaN
represented a missing value. But in Dune, I think that's something we hardly ever want. If a vector produced by a computation contains NaN
somewhere, it typically means something went very wrong and you want DenseVector::infinity_norm()
to tell you that just like DenseVector::two_norm()
would. So the return value should be NaN
if any entry is NaN
.
I've raised this issue before. And we did something about it: Clearly, std::fmax
is not a solution; std::max
is in a way even worse because what it returns depends on the order that you pass your arguments in (i.e. std::fmax(3, std::nan("")) == 3 == std::fmax(std::nan(""), 3)
while no such thing holds for std::max()
). One could of course manually check each entry via std::is_nan()
but that is expensive. Too expensive, even. So in flyspray/FS#1147 (closed) it was decided that we go with a construct that looks a little like this:
// compute the actual return value ret (turned one loop into two loops here for better readability)
...
real_type isNaN = 1;
for (auto const &x : *this) {
real_type const a = abs(x);
isNaN += a;
}
isNaN /= isNaN;
return ret * isNaN;
mostly for reasons of performance. The logic being: If any x
is NaN
, then a = abs(x)
is NaN
, thus isNaN += a
makes isNaN
into NaN
. Finally, isNaN /= isNaN
will make isNaN
either 1
or NaN
, so that diving ret
by it will either leave ret
unchanged or turn it into NaN
.
Now here's the problem: While this chain of ifs is true, they're not if and only ifs. In particular, what if x
is std::numeric_limits<double>::infinity()
(I'll just write inf
for short). If x
is inf
, then a
is inf, then (if isNaN
was a finite number) isNaN += a
will make isNaN
into inf
and then isNaN /= isNaN
turns isNaN
into... NaN
.
TL;DR: A vector that contains inf
will currently produce an infinity_norm
that is NaN
rather than inf
(actually, it returns -NaN
rather than NaN
but I'm not sure why), while two_norm()
will return inf
.
So I'm afraid we should revisit @christi's list of proposed implementations from here and pick a slightly slower but safer one.