Merge branch 'feature/FS1147-norm-of-vectors-with-nan' into 'master'

flyspray/FS#1147: Compute norm of vectors with NaN entries properly

This is an attempt to solve the issue outlined in flyspray/FS#1147. I've yet to fix `DenseMatrix::infinity_norm` the same way.

See merge request !15

dune/common/densematrix.hh

@@ -567,31 +567,69 @@ namespace Dune
     }
 
     //! infinity norm (row sum norm, how to generalize for blocks?)
-    typename FieldTraits<value_type>::real_type infinity_norm () const
-    {
-      if (size() == 0)
-        return 0.0;
-
-      ConstIterator it = begin();
-      typename remove_const< typename FieldTraits<value_type>::real_type >::type max = (*it).one_norm();
-      for (it = it + 1; it != end(); ++it)
-        max = std::max(max, (*it).one_norm());
-
-      return max;
+    template <typename vt = value_type,
+              typename std::enable_if<!has_nan<vt>::value, int>::type = 0>
+    typename FieldTraits<vt>::real_type infinity_norm() const {
+      using real_type = typename FieldTraits<vt>::real_type;
+      using std::max;
+
+      real_type norm = 0;
+      for (auto const &x : *this) {
+        real_type const a = x.one_norm();
+        norm = max(a, norm);
+      }
+      return norm;
     }
 
     //! simplified infinity norm (uses Manhattan norm for complex values)
-    typename FieldTraits<value_type>::real_type infinity_norm_real () const
-    {
-      if (size() == 0)
-        return 0.0;
+    template <typename vt = value_type,
+              typename std::enable_if<!has_nan<vt>::value, int>::type = 0>
+    typename FieldTraits<vt>::real_type infinity_norm_real() const {
+      using real_type = typename FieldTraits<vt>::real_type;
+      using std::max;
+
+      real_type norm = 0;
+      for (auto const &x : *this) {
+        real_type const a = x.one_norm_real();
+        norm = max(a, norm);
+      }
+      return norm;
+    }
 
-      ConstIterator it = begin();
-      typename remove_const< typename FieldTraits<value_type>::real_type >::type max = (*it).one_norm_real();
-      for (it = it + 1; it != end(); ++it)
-        max = std::max(max, (*it).one_norm_real());
+    //! infinity norm (row sum norm, how to generalize for blocks?)
+    template <typename vt = value_type,
+              typename std::enable_if<has_nan<vt>::value, int>::type = 0>
+    typename FieldTraits<vt>::real_type infinity_norm() const {
+      using real_type = typename FieldTraits<vt>::real_type;
+      using std::max;
+
+      real_type norm = 0;
+      real_type isNaN = 1;
+      for (auto const &x : *this) {
+        real_type const a = x.one_norm();
+        norm = max(a, norm);
+        isNaN += a;
+      }
+      isNaN /= isNaN;
+      return norm * isNaN;
+    }
 
-      return max;
+    //! simplified infinity norm (uses Manhattan norm for complex values)
+    template <typename vt = value_type,
+              typename std::enable_if<has_nan<vt>::value, int>::type = 0>
+    typename FieldTraits<vt>::real_type infinity_norm_real() const {
+      using real_type = typename FieldTraits<vt>::real_type;
+      using std::max;
+
+      real_type norm = 0;
+      real_type isNaN = 1;
+      for (auto const &x : *this) {
+        real_type const a = x.one_norm_real();
+        norm = max(a, norm);
+        isNaN += a;
+      }
+      isNaN /= isNaN;
+      return norm * isNaN;
     }
 
     //===== solve

dune/common/densevector.hh

@@ -604,37 +604,71 @@ namespace Dune {
     }
 
     //! infinity norm (maximum of absolute values of entries)
-    typename FieldTraits<value_type>::real_type infinity_norm () const
-    {
+    template <typename vt = value_type,
+              typename std::enable_if<!has_nan<vt>::value, int>::type = 0>
+    typename FieldTraits<vt>::real_type infinity_norm() const {
+      using real_type = typename FieldTraits<vt>::real_type;
       using std::abs;
       using std::max;
-      typedef typename FieldTraits<value_type>::real_type real_type;
 
-      if (size() == 0)
-        return 0.0;
+      real_type norm = 0;
+      for (auto const &x : *this) {
+        real_type const a = abs(x);
+        norm = max(a, norm);
+      }
+      return norm;
+    }
 
-      ConstIterator it = begin();
-      real_type max_val = abs(*it);
-      for (it = it + 1; it != end(); ++it)
-        max_val = max(max_val, real_type(abs(*it)));
+    //! simplified infinity norm (uses Manhattan norm for complex values)
+    template <typename vt = value_type,
+              typename std::enable_if<!has_nan<vt>::value, int>::type = 0>
+    typename FieldTraits<vt>::real_type infinity_norm_real() const {
+      using real_type = typename FieldTraits<vt>::real_type;
+      using std::max;
 
-      return max_val;
+      real_type norm = 0;
+      for (auto const &x : *this) {
+        real_type const a = fvmeta::absreal(x);
+        norm = max(a, norm);
+      }
+      return norm;
     }
 
-    //! simplified infinity norm (uses Manhattan norm for complex values)
-    typename FieldTraits<value_type>::real_type infinity_norm_real () const
-    {
+    //! infinity norm (maximum of absolute values of entries)
+    template <typename vt = value_type,
+              typename std::enable_if<has_nan<vt>::value, int>::type = 0>
+    typename FieldTraits<vt>::real_type infinity_norm() const {
+      using real_type = typename FieldTraits<vt>::real_type;
+      using std::abs;
       using std::max;
 
-      if (size() == 0)
-        return 0.0;
+      real_type norm = 0;
+      real_type isNaN = 1;
+      for (auto const &x : *this) {
+        real_type const a = abs(x);
+        norm = max(a, norm);
+        isNaN += a;
+      }
+      isNaN /= isNaN;
+      return norm * isNaN;
+    }
 
-      ConstIterator it = begin();
-      typename FieldTraits<value_type>::real_type max_val = fvmeta::absreal(*it);
-      for (it = it + 1; it != end(); ++it)
-        max_val = max(max_val, fvmeta::absreal(*it));
+    //! simplified infinity norm (uses Manhattan norm for complex values)
+    template <typename vt = value_type,
+              typename std::enable_if<has_nan<vt>::value, int>::type = 0>
+    typename FieldTraits<vt>::real_type infinity_norm_real() const {
+      using real_type = typename FieldTraits<vt>::real_type;
+      using std::max;
 
-      return max_val;
+      real_type norm = 0;
+      real_type isNaN = 1;
+      for (auto const &x : *this) {
+        real_type const a = fvmeta::absreal(x);
+        norm = max(a, norm);
+        isNaN += a;
+      }
+      isNaN /= isNaN;
+      return norm * isNaN;
     }
 
     //===== sizes

dune/common/test/fmatrixtest.cc

@@ -502,27 +502,64 @@ void test_invert ()
   A.invert();
 }
 
-// Make sure that a matrix with only NaN entries has norm NaN.
-// Prior to r6819, the infinity_norm would be zero; see also FS #1147.
+template <class M>
+void checkNormNAN(M const &v, int line) {
+  if (!std::isnan(v.frobenius_norm())) {
+    std::cerr << "error: norm not NaN: frobenius_norm() on line "
+              << line << " (type: " << Dune::className(v[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]) << ")"
+              << std::endl;
+    std::exit(-1);
+  }
+}
+
+// Make sure that matrices with NaN entries have norm NaN.
+// See also bug flyspray/FS#1147
+template <typename T>
 void
-test_nan()
+test_nan(T const &mynan)
 {
-  double mynan = 0.0/0.0;
-
-  Dune::FieldMatrix<double, 2, 2> m2(mynan);
-  assert(std::isnan(m2.infinity_norm()));
-  assert(std::isnan(m2.frobenius_norm()));
-  assert(std::isnan(m2.frobenius_norm2()));
-
-  Dune::FieldMatrix<double, 0, 2> m02(mynan);
-  assert(0.0 == m02.infinity_norm());
-  assert(0.0 == m02.frobenius_norm());
-  assert(0.0 == m02.frobenius_norm2());
-
-  Dune::FieldMatrix<double, 2, 0> m20(mynan);
-  assert(0.0 == m20.infinity_norm());
-  assert(0.0 == m20.frobenius_norm());
-  assert(0.0 == m20.frobenius_norm2());
+  T const n(0);
+  {
+    Dune::FieldMatrix<T, 2, 2> m = {
+      { mynan, mynan },
+      { mynan, mynan }
+    };
+    checkNormNAN(m, __LINE__);
+  }
+  {
+    Dune::FieldMatrix<T, 2, 2> m = {
+      { mynan, n },
+      { n, n }
+    };
+    checkNormNAN(m, __LINE__);
+  }
+  {
+    Dune::FieldMatrix<T, 2, 2> m = {
+      { n, mynan },
+      { n, n }
+    };
+    checkNormNAN(m, __LINE__);
+  }
+  {
+    Dune::FieldMatrix<T, 2, 2> m = {
+      { n, n },
+      { mynan, n }
+    };
+    checkNormNAN(m, __LINE__);
+  }
+  {
+    Dune::FieldMatrix<T, 2, 2> m = {
+      { n, n },
+      { n, mynan }
+    };
+    checkNormNAN(m, __LINE__);
+  }
 }
 
 // The computation of infinity_norm_real() was flawed from r6819 on
@@ -568,7 +605,14 @@ void test_initialisation()
 int main()
 {
   try {
-    test_nan();
+    {
+      double nan = std::nan("");
+      test_nan(nan);
+    }
+    {
+      std::complex<double> nan( std::nan(""), 17 );
+      test_nan(nan);
+    }
     test_infinity_norms();
     test_initialisation();
 

dune/common/test/fvectortest.cc

@@ -425,24 +425,46 @@ public:
   }
 };
 
-// Make sure that a vector with only NaN entries has norm NaN.
-// Prior to r6914, the infinity_norm would be zero; see also FS #1147.
+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]) << ")"
+              << 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]) << ")"
+              << 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]) << ")"
+              << 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()
+test_nan(T const &mynan)
 {
-  double mynan = 0.0/0.0;
-
-  Dune::FieldVector<double, 2> v2(mynan);
-  assert(std::isnan(v2.infinity_norm()));
-  assert(std::isnan(v2.one_norm()));
-  assert(std::isnan(v2.two_norm()));
-  assert(std::isnan(v2.two_norm2()));
-
-  Dune::FieldVector<double, 0> v0(mynan);
-  assert(0.0 == v0.infinity_norm());
-  assert(0.0 == v0.one_norm());
-  assert(0.0 == v0.two_norm());
-  assert(0.0 == v0.two_norm2());
+  {
+    Dune::FieldVector<T, 2> v = { mynan, mynan };
+    checkNormNAN(v, __LINE__);
+  }
+  {
+    Dune::FieldVector<T, 2> v = { mynan, 0 };
+    checkNormNAN(v, __LINE__);
+  }
+  {
+    Dune::FieldVector<T, 2> v = { 0, mynan };
+    checkNormNAN(v, __LINE__);
+  }
 }
 
 void
@@ -486,7 +508,14 @@ int main()
     DotProductTest<ft,3>();
 #endif // HAVE_GMP
 
-    test_nan();
+    {
+      double nan = std::nan("");
+      test_nan(nan);
+    }
+    {
+      std::complex<double> nan( std::nan(""), 17 );
+      test_nan(nan);
+    }
     test_infinity_norms();
     test_initialisation();
   } catch (Dune::Exception& e) {

dune/common/typetraits.hh

@@ -405,6 +405,15 @@ namespace Dune
     static const bool value = true;
   };
 
+  template <typename T>
+  struct has_nan
+      : public std::integral_constant<bool, std::is_floating_point<T>::value> {
+  };
+
+  template <typename T>
+  struct has_nan<std::complex<T>>
+      : public std::integral_constant<bool, std::is_floating_point<T>::value> {
+  };
 
 #if defined(DOXYGEN) or HAVE_IS_INDEXABLE_SUPPORT