fixed min/max bug in FieldMatVec

src/amdis/common/FieldMatVec.hpp

 #pragma once
 
 #include <algorithm>
+#include <limits>
 
 #include <dune/common/diagonalmatrix.hh>
 #include <dune/common/fmatrix.hh>
@@ -96,21 +97,19 @@ namespace AMDiS
   namespace Impl
   {
     template <class T, int N, class Operation>
-    T accumulate(FieldVector<T, N> const& x, Operation op)
+    T accumulate(FieldVector<T, N> const& x, T init, Operation op)
     {
-      T result = 0;
       for (int i = 0; i < N; ++i)
-        result = op(result, x[i]);
-      return result;
+        init = op(init, x[i]);
+      return init;
     }
 
     template <class T, int N, class Operation>
-    T accumulate(FieldMatrix<T, 1, N> const& x, Operation op)
+    T accumulate(FieldMatrix<T, 1, N> const& x, T init, Operation op)
     {
-      T result = 0;
       for (int i = 0; i < N; ++i)
-        result = op(result, x[0][i]);
-      return result;
+        init = op(init, x[0][i]);
+      return init;
     }
 
   } // end namespace Impl
@@ -119,13 +118,13 @@ namespace AMDiS
   template <class T, int N>
   T sum(FieldVector<T, N> const& x)
   {
-    return Impl::accumulate(x, Operation::Plus{});
+    return Impl::accumulate(x, T(0), Operation::Plus{});
   }
 
   template <class T, int N>
   T sum(FieldMatrix<T, 1, N> const& x)
   {
-    return Impl::accumulate(x, Operation::Plus{});
+    return Impl::accumulate(x, T(0), Operation::Plus{});
   }
 
 
@@ -134,66 +133,66 @@ namespace AMDiS
   auto unary_dot(FieldVector<T, N> const& x)
   {
     auto op = [](auto const& a, auto const& b) { return a + Math::sqr(std::abs(b)); };
-    return Impl::accumulate(x, op);
+    return Impl::accumulate(x, T(0), op);
   }
 
   template <class T, int N>
   auto unary_dot(FieldMatrix<T, 1, N> const& x)
   {
     auto op = [](auto const& a, auto const& b) { return a + Math::sqr(std::abs(b)); };
-    return Impl::accumulate(x, op);
+    return Impl::accumulate(x, T(0), op);
   }
 
   /// Maximum over all vector entries
   template <class T, int N>
   auto max(FieldVector<T, N> const& x)
   {
-    return Impl::accumulate(x, Operation::Max{});
+    return Impl::accumulate(x, std::numeric_limits<T>::lowest(), Operation::Max{});
   }
 
   template <class T, int N>
   auto max(FieldMatrix<T, 1, N> const& x)
   {
-    return Impl::accumulate(x, Operation::Max{});
+    return Impl::accumulate(x, std::numeric_limits<T>::lowest(), Operation::Max{});
   }
 
   /// Minimum over all vector entries
   template <class T, int N>
   auto min(FieldVector<T, N> const& x)
   {
-    return Impl::accumulate(x, Operation::Min{});
+    return Impl::accumulate(x, std::numeric_limits<T>::max(), Operation::Min{});
   }
 
   template <class T, int N>
   auto min(FieldMatrix<T, 1, N> const& x)
   {
-    return Impl::accumulate(x, Operation::Min{});
+    return Impl::accumulate(x, std::numeric_limits<T>::max(), Operation::Min{});
   }
 
   /// Maximum of the absolute values of vector entries
   template <class T, int N>
   auto abs_max(FieldVector<T, N> const& x)
   {
-    return Impl::accumulate(x, Operation::AbsMax{});
+    return Impl::accumulate(x, T(0), Operation::AbsMax{});
   }
 
   template <class T, int N>
   auto abs_max(FieldMatrix<T, 1, N> const& x)
   {
-    return Impl::accumulate(x, Operation::AbsMax{});
+    return Impl::accumulate(x, T(0), Operation::AbsMax{});
   }
 
   /// Minimum of the absolute values of vector entries
   template <class T, int N>
   auto abs_min(FieldVector<T, N> const& x)
   {
-    return Impl::accumulate(x, Operation::AbsMin{});
+    return Impl::accumulate(x, std::numeric_limits<T>::max(), Operation::AbsMin{});
   }
 
   template <class T, int N>
   auto abs_min(FieldMatrix<T, 1, N> const& x)
   {
-    return Impl::accumulate(x, Operation::AbsMin{});
+    return Impl::accumulate(x, std::numeric_limits<T>::max(), Operation::AbsMin{});
   }
 
   // ----------------------------------------------------------------------------
@@ -205,14 +204,14 @@ namespace AMDiS
   auto one_norm(FieldVector<T, N> const& x)
   {
     auto op = [](auto const& a, auto const& b) { return a + std::abs(b); };
-    return Impl::accumulate(x, op);
+    return Impl::accumulate(x, T(0), op);
   }
 
   template <class T, int N>
   auto one_norm(FieldMatrix<T, 1, N> const& x)
   {
     auto op = [](auto const& a, auto const& b) { return a + std::abs(b); };
-    return Impl::accumulate(x, op);
+    return Impl::accumulate(x, T(0), op);
   }
 
   /** \ingroup vector_norms
@@ -237,14 +236,14 @@ namespace AMDiS
   auto p_norm(FieldVector<T, N> const& x)
   {
     auto op = [](auto const& a, auto const& b) { return a + Math::pow<p>(std::abs(b)); };
-    return std::pow( Impl::accumulate(x, op), 1.0/p );
+    return std::pow( Impl::accumulate(x, T(0), op), 1.0/p );
   }
 
   template <int p, class T, int N>
   auto p_norm(FieldMatrix<T, 1, N> const& x)
   {
     auto op = [](auto const& a, auto const& b) { return a + Math::pow<p>(std::abs(b)); };
-    return std::pow( Impl::accumulate(x, op), 1.0/p );
+    return std::pow( Impl::accumulate(x, T(0), op), 1.0/p );
   }
 
   /** \ingroup vector_norms

test/CMakeLists.txt

 
+dune_add_test(SOURCES ClonablePtrTest.cpp
+  LINK_LIBRARIES amdis)
+  
 dune_add_test(SOURCES ConceptsTest.cpp
   LINK_LIBRARIES amdis)
 
@@ -29,3 +32,6 @@ dune_add_test(SOURCES StringTest.cpp
 
 dune_add_test(SOURCES TreeDataTest.cpp
   LINK_LIBRARIES amdis)
+
+dune_add_test(SOURCES TupleUtilityTest.cpp
+  LINK_LIBRARIES amdis)

test/FieldMatVecTest.cpp

@@ -53,6 +53,18 @@ void test1()
   AMDIS_TEST_EQ( max(d), -1.0 );
   AMDIS_TEST_EQ( abs_min(c), 3.0 );
   AMDIS_TEST_EQ( abs_max(c), 5.0 );
+
+  using V3 = FieldVector<int, 3>;
+  V3 f{2, 3, 4};
+  V3 g{3, -5, 4};
+  V3 h{-1, -2, -3};
+  AMDIS_TEST_EQ( sum(g),  2 );
+  AMDIS_TEST_EQ( min(f),  2 );
+  AMDIS_TEST_EQ( min(g), -5 );
+  AMDIS_TEST_EQ( max(g),  4 );
+  AMDIS_TEST_EQ( max(h), -1 );
+  AMDIS_TEST_EQ( abs_min(g), 3 );
+  AMDIS_TEST_EQ( abs_max(g), 5 );
 }
 
 // -----------------------------------------------------------------------------