diff --git a/dune/istl/bvector.hh b/dune/istl/bvector.hh index fa2251a93b9175a1dfeff0dfc26cca3438a51996..e22a4fd7f572329bb0603647d124648f753de318 100644 --- a/dune/istl/bvector.hh +++ b/dune/istl/bvector.hh @@ -196,19 +196,69 @@ namespace Dune { } //! infinity norm (maximum of absolute values of entries) - typename FieldTraits<field_type>::real_type infinity_norm () const - { - typename FieldTraits<field_type>::real_type max=0; - for (size_type i=0; i<this->n; ++i) max = std::max(max,(*this)[i].infinity_norm()); - return max; + template <typename ft = field_type, + typename std::enable_if<!has_nan<ft>::value, int>::type = 0> + typename FieldTraits<ft>::real_type infinity_norm() const { + using real_type = typename FieldTraits<ft>::real_type; + using std::max; + + real_type norm = 0; + for (auto const &x : *this) { + real_type const a = x.infinity_norm(); + norm = max(a, norm); + } + return norm; } //! simplified infinity norm (uses Manhattan norm for complex values) - typename FieldTraits<field_type>::real_type infinity_norm_real () const - { - typename FieldTraits<field_type>::real_type max=0; - for (size_type i=0; i<this->n; ++i) max = std::max(max,(*this)[i].infinity_norm_real()); - return max; + template <typename ft = field_type, + typename std::enable_if<!has_nan<ft>::value, int>::type = 0> + typename FieldTraits<ft>::real_type infinity_norm_real() const { + using real_type = typename FieldTraits<ft>::real_type; + using std::max; + + real_type norm = 0; + for (auto const &x : *this) { + real_type const a = x.infinity_norm_real(); + norm = max(a, norm); + } + return norm; + } + + //! infinity norm (maximum of absolute values of entries) + template <typename ft = field_type, + typename std::enable_if<has_nan<ft>::value, int>::type = 0> + typename FieldTraits<ft>::real_type infinity_norm() const { + using real_type = typename FieldTraits<ft>::real_type; + using std::max; + + real_type norm = 0; + real_type isNaN = 1; + for (auto const &x : *this) { + real_type const a = x.infinity_norm(); + norm = max(a, norm); + isNaN += a; + } + isNaN /= isNaN; + return norm * isNaN; + } + + //! simplified infinity norm (uses Manhattan norm for complex values) + template <typename ft = field_type, + typename std::enable_if<has_nan<ft>::value, int>::type = 0> + typename FieldTraits<ft>::real_type infinity_norm_real() const { + using real_type = typename FieldTraits<ft>::real_type; + using std::max; + + real_type norm = 0; + real_type isNaN = 1; + for (auto const &x : *this) { + real_type const a = x.infinity_norm_real(); + norm = max(a, norm); + isNaN += a; + } + isNaN /= isNaN; + return norm * isNaN; } //===== sizes diff --git a/dune/istl/test/bvectortest.cc b/dune/istl/test/bvectortest.cc index ca7fdce2589246b0ab1f4a95b4946945c9a95e8b..5e42a99956c8d590bdd6197650ed3adb7d214dd6 100644 --- a/dune/istl/test/bvectortest.cc +++ b/dune/istl/test/bvectortest.cc @@ -5,6 +5,7 @@ #include <dune/common/fvector.hh> #include <dune/common/poolallocator.hh> #include <dune/common/debugallocator.hh> +#include <dune/common/classname.hh> template<typename T, int BS> void assign(Dune::FieldVector<T,BS>& b, const T& i) @@ -120,6 +121,74 @@ void testCapacity() vec1.reserve(0, false); } +template <class V> +void checkNormNAN(V const &v, int line) { + if (!std::isnan(v.one_norm())) { + std::cerr << "error: norm not NaN: one_norm() on line " + << line << " (type: " << Dune::className(v[0][0]) << ")" + << std::endl; + std::exit(-1); + } + if (!std::isnan(v.two_norm())) { + std::cerr << "error: norm not NaN: two_norm() on line " + << line << " (type: " << Dune::className(v[0][0]) << ")" + << std::endl; + std::exit(-1); + } + if (!std::isnan(v.infinity_norm())) { + std::cerr << "error: norm not NaN: infinity_norm() on line " + << line << " (type: " << Dune::className(v[0][0]) << ")" + << std::endl; + std::exit(-1); + } +} + +// Make sure that vectors with NaN entries have norm NaN. +// See also bug flyspray/FS#1147 +template <typename T> +void +test_nan(T const &mynan) +{ + using FV = Dune::FieldVector<T,2>; + using V = Dune::BlockVector<FV>; + T n(0); + { + V v = { + { mynan, n }, + { n, n } + }; + checkNormNAN(v, __LINE__); + } + { + V v = { + { n, mynan }, + { n, n } + }; + checkNormNAN(v, __LINE__); + } + { + V v = { + { n, n }, + { mynan, n } + }; + checkNormNAN(v, __LINE__); + } + { + V v = { + { n, n }, + { n, mynan } + }; + checkNormNAN(v, __LINE__); + } + { + V v = { + { mynan, mynan }, + { mynan, mynan } + }; + checkNormNAN(v, __LINE__); + } +} + int main() { typedef std::complex<double> value_type; @@ -141,6 +210,15 @@ int main() assert(fromInitializerList[1] == value_type(1)); assert(fromInitializerList[2] == value_type(2)); + { + double nan = std::nan(""); + test_nan(nan); + } + { + std::complex<double> nan( std::nan(""), 17 ); + test_nan(nan); + } + int ret = 0; ret += testVector<1>();