## 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 **if**s is true, they're not **if and only if**s. 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.