Implemented float comparison with templates

CMakeLists_files.cmake

@@ -22,7 +22,6 @@
 
 list (APPEND MAIN_SOURCE_FILES
       opm/common/data/SimulationDataContainer.cpp
-      opm/common/util/numeric/cmp.cpp
 )
 
 list (APPEND TEST_SOURCE_FILES

opm/common/data/SimulationDataContainer.cpp

@@ -119,7 +119,7 @@ namespace Opm {
 
             if (other.hasCellData( key )) {
                 const auto& other_data = other.getCellData( key );
-                if (!cmp::double_vector_equal( data , other_data ))
+                if (!cmp::vector_equal<double>( data , other_data ))
                     return false;
             } else
                 return false;
@@ -131,7 +131,7 @@ namespace Opm {
 
             if (other.hasFaceData( key )) {
                 const auto& other_data = other.getFaceData( key );
-                if (!cmp::double_vector_equal( data , other_data ))
+                if (!cmp::vector_equal<double>( data , other_data ))
                     return false;
             } else
                 return false;

opm/common/util/numeric/cmp.cpp

@@ -1,91 +0,0 @@
-/*
-  Copyright 2016 Statoil ASA.
-
-  This file is part of the Open Porous Media project (OPM).
-
-  OPM is free software: you can redistribute it and/or modify
-  it under the terms of the GNU General Public License as published by
-  the Free Software Foundation, either version 3 of the License, or
-  (at your option) any later version.
-
-  OPM is distributed in the hope that it will be useful,
-  but WITHOUT ANY WARRANTY; without even the implied warranty of
-  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-  GNU General Public License for more details.
-
-  You should have received a copy of the GNU General Public License
-  along with OPM.  If not, see <http://www.gnu.org/licenses/>.
-*/
-
-#include <math.h>
-#include <string.h>
-
-#include <algorithm>
-#include <cstddef>
-#include <vector>
-
-#include <opm/common/util/numeric/cmp.hpp>
-
-
-namespace Opm {
-    namespace cmp {
-
-        bool double_equal(double value1 , double value2 , double abs_eps , double rel_eps) {
-            bool equal = true;
-            double diff = fabs(value1 - value2);
-            if (diff > abs_eps) {
-                double scale = std::max(fabs(value1), fabs(value2));
-
-                if (diff > scale * rel_eps) {
-                    equal = false;
-                }
-            }
-            return equal;
-        }
-
-
-        bool double_equal(double value1 , double value2) {
-            return double_equal( value1 , value2 , default_abs_epsilon , default_rel_epsilon );
-        }
-
-        /*****************************************************************/
-
-        bool double_vector_equal(const std::vector<double>& v1, const std::vector<double>& v2, double abs_eps, double rel_eps) {
-            if (v1.size() != v2.size()) {
-                return false;
-            }
-
-            for (size_t i = 0; i < v1.size(); i++) {
-                if (!double_equal( v1[i] , v2[i]))
-                    return false;
-            }
-
-            return true;
-        }
-
-        bool double_vector_equal(const std::vector<double>& v1, const std::vector<double>& v2) {
-            return double_vector_equal( v1 , v2 , default_abs_epsilon , default_rel_epsilon);
-        }
-
-
-        bool double_ptr_equal(const double * p1, const double *p2, size_t num_elements, double abs_eps , double rel_eps) {
-            if (memcmp(p1 , p2 , num_elements * sizeof * p1) == 0)
-                return true;
-            else {
-                size_t index;
-                for (index = 0; index < num_elements; index++) {
-                    if (!double_equal( p1[index] , p2[index] , abs_eps , rel_eps)) {
-                        return false;
-                    }
-                }
-            }
-            return true;
-        }
-
-
-        bool double_ptr_equal(const double * p1, const double *p2, size_t num_elements) {
-            return double_ptr_equal( p1 , p2  ,num_elements , default_abs_epsilon , default_rel_epsilon );
-        }
-
-    }
-}

opm/common/util/numeric/cmp.hpp

@@ -22,6 +22,8 @@
 
 #include <cstddef>
 #include <vector>
+#include <type_traits>
+#include <math.h>
 
 namespace Opm {
 
@@ -29,38 +31,101 @@ namespace Opm {
     /// approximate comparison of double values based on absolute
     /// and relative difference. There are three functions:
     ///
-    ///   double_equal() : Compare two double values.
+    ///   scalar_equal<T>() : Compare two <T> values.
     ///
-    ///   double_ptr_equal(): This compares all the element in the
+    ///   ptr_equal<T>(): This compares all the element in the
-    ///      two double * pointers.
+    ///      two T * pointers.
     ///
-    ///   double_vector_equal(): This compares all the elements in
+    ///   vector_equal<T>(): This compares all the elements in
-    ///      two std::vector<double> instances.
+    ///      two std::vector<T> instances.
     ///
-    /// For both double_vector_equal() and double_ptr_equal() the
+    /// For both vector_equal<T>() and ptr_equal<T>() the
-    /// actual comparison is based on the double_equal()
+    /// actual comparison is based on the scalar_equal<T>()
     /// function. All functions exist as two overloads, one which
     /// takes explicit input values for the absolute and relative
-    /// epsilon, and one which uses default values.
+    /// epsilon, and one which uses default values.  
+    ///
+    /// The comparison functions are implemented as templates, with
+    /// the following caveats:
+    ///
+    ///  1. The static_assert() in scalar_equal<T> ensures that only
+    ///     floating point types can be used.
+    ///
+    ///  2. The default epsilon values are of type double -
+    ///     irrespective of the type of data being compared.
     ///
     /// For more details of floating point comparison please consult
     /// this reference:
     ///
     ///    https://randomascii.wordpress.com/2012/02/25/comparing-floating-point-numbers-2012-edition/
+
     namespace cmp {
 
         const double default_abs_epsilon = 1e-8;
         const double default_rel_epsilon = 1e-5;
 
+        template<typename T>
+        bool scalar_equal(T value1, T value2, T abs_eps , T rel_eps) {
+            static_assert(std::is_floating_point<T>::value, "Function scalar_equal() A can only be instantiated with floating point types");
 
-        bool double_equal(double value1, double value2, double abs_eps , double rel_eps);
+	    bool equal = true;
-        bool double_equal(double value1, double value2);
+            T diff = fabs(value1 - value2);
+            if (diff > abs_eps) {
+                T scale = std::max(fabs(value1), fabs(value2));
 
-        bool double_vector_equal(const std::vector<double>& v1, const std::vector<double>& v2, double abs_eps, double rel_eps);
+                if (diff > scale * rel_eps) {
-        bool double_vector_equal(const std::vector<double>& v1, const std::vector<double>& v2);
+                    equal = false;
+		}
+	    }
+            return equal;
+	}
 
-        bool double_ptr_equal(const double* p1, const double* p2, size_t num_elements, double abs_eps, double rel_eps);
+
-        bool double_ptr_equal(const double* p1, const double* p2, size_t num_elements);
+        template<typename T>
+        bool scalar_equal(T value1, T value2) {
+  	    return scalar_equal<T>( value1 , value2 , default_abs_epsilon , default_rel_epsilon );
+	}
+
+        template<typename T>
+        bool vector_equal(const std::vector<T>& v1, const std::vector<T>& v2, T abs_eps, T rel_eps) {
+          if (v1.size() != v2.size()) {
+                return false;
+            }
+
+            for (size_t i = 0; i < v1.size(); i++) {
+                if (!scalar_equal<T>( v1[i] , v2[i]))
+                    return false;
+            }
+
+            return true;
+        }
+
+        template<typename T>
+        bool vector_equal(const std::vector<T>& v1, const std::vector<T>& v2) {
+            return vector_equal<T>(v1, v2, default_abs_epsilon, default_rel_epsilon);
+        }
+
+
+        template<typename T>
+        bool ptr_equal(const T* p1, const T* p2, size_t num_elements, T abs_eps, T rel_eps) {
+            if (memcmp(p1 , p2 , num_elements * sizeof * p1) == 0)
+                return true;
+            else {
+                size_t index;
+                for (index = 0; index < num_elements; index++) {
+                    if (!scalar_equal<T>( p1[index] , p2[index] , abs_eps , rel_eps)) {
+                        return false;
+                    }
+                }
+            }
+            return true;
+        }
+
+
+        template<typename T>
+        bool ptr_equal(const T* p1, const T* p2, size_t num_elements) {
+             return ptr_equal<T>(p1, p2, num_elements , default_abs_epsilon, default_rel_epsilon);
+        }
     }
 }
 

tests/test_cmp.cpp

@@ -38,10 +38,10 @@ BOOST_AUTO_TEST_CASE(TestSCalarcmp) {
     const double abs_epsilon = cmp::default_abs_epsilon;
     const double rel_epsilon = cmp::default_rel_epsilon;
 
-    BOOST_CHECK( cmp::double_equal(1,1));
+    BOOST_CHECK( cmp::scalar_equal<double>(1,1));
-    BOOST_CHECK_EQUAL( false , cmp::double_equal(1,0));
+    BOOST_CHECK_EQUAL( false , cmp::scalar_equal<double>(1,0));
-    BOOST_CHECK_EQUAL( false , cmp::double_equal(0,1));
+    BOOST_CHECK_EQUAL( false , cmp::scalar_equal<double>(0,1));
-    BOOST_CHECK_EQUAL( false , cmp::double_equal(-1,1));
+    BOOST_CHECK_EQUAL( false , cmp::scalar_equal<double>(-1,1));
 
 
 
@@ -50,60 +50,74 @@ BOOST_AUTO_TEST_CASE(TestSCalarcmp) {
     {
         v1 = 0.0;
         v2 = 0.0;
-        BOOST_CHECK( cmp::double_equal( v1 , v2));
+        BOOST_CHECK( cmp::scalar_equal<double>( v1 , v2));
 
         v1 = 1e-12;
         v2 = v1 + 0.5*abs_epsilon;
-        BOOST_CHECK( cmp::double_equal( v1 , v2));
+        BOOST_CHECK( cmp::scalar_equal<double>( v1 , v2));
 
         v1 = 7.0;
         v2 = 7.0;
-        BOOST_CHECK( cmp::double_equal( v1 , v2));
+        BOOST_CHECK( cmp::scalar_equal<double>( v1 , v2));
 
         v1 = -7.0;
         v2 = -7.0;
-        BOOST_CHECK( cmp::double_equal( v1 , v2));
+        BOOST_CHECK( cmp::scalar_equal<double>( v1 , v2));
 
         v1 = 0;
         v2 = 0.5 * abs_epsilon;
-        BOOST_CHECK( cmp::double_equal( v1 , v2));
+        BOOST_CHECK( cmp::scalar_equal<double>( v1 , v2));
 
 
         v1 = 1e7;
         v2 = 1e7 + 2*abs_epsilon;
-        BOOST_CHECK( cmp::double_equal( v1 , v2 ));
+        BOOST_CHECK( cmp::scalar_equal<double>( v1 , v2 ));
 
         v1 = 1e7*(1 - abs_epsilon);
         v2 = 1e7*(1 + rel_epsilon);
-        BOOST_CHECK( !cmp::double_equal( v1 , v2 ));
+        BOOST_CHECK( !cmp::scalar_equal<double>( v1 , v2 ));
 
         v1 = 1e7*(1 + abs_epsilon);
         v2 = 1e7*(1 + rel_epsilon);
-        BOOST_CHECK( cmp::double_equal( v1 , v2 ));
+        BOOST_CHECK( cmp::scalar_equal<double>( v1 , v2 ));
     }
 
     /* Should be different: */
     {
         v1 = 0;
         v2 = 1.5 * abs_epsilon;
-        BOOST_CHECK( !cmp::double_equal( v1 , v2 ));
+        BOOST_CHECK( !cmp::scalar_equal<double>( v1 , v2 ));
 
         v1 = 1e-8;
         v2 = v1 + 1.5*abs_epsilon;
-        BOOST_CHECK( !cmp::double_equal( v1 , v2 ));
+        BOOST_CHECK( !cmp::scalar_equal<double>( v1 , v2 ));
 
         v1 = 1;
         v2 = v1*(1 + 2*rel_epsilon + abs_epsilon);
-        BOOST_CHECK( !cmp::double_equal( v1 , v2 ));
+        BOOST_CHECK( !cmp::scalar_equal<double>( v1 , v2 ));
 
         v1 = 10;
         v2 = v1*(1 + 2*rel_epsilon + abs_epsilon);
-        BOOST_CHECK( !cmp::double_equal( v1 , v2 ));
+        BOOST_CHECK( !cmp::scalar_equal<double>( v1 , v2 ));
 
         v1 = 1e7;
         v2 = 1e7*(1 + 2*rel_epsilon + abs_epsilon);
-        BOOST_CHECK( !cmp::double_equal( v1 , v2 ));
+        BOOST_CHECK( !cmp::scalar_equal<double>( v1 , v2 ));
     }
 }
 
+/* Ensure that float instantiation works. */
+BOOST_AUTO_TEST_CASE(TestFloatcmp) {
+    std::vector<float> v1;
+    std::vector<float> v2;
+    for (size_t i =0; i < 10; i++) {
+        v1.push_back( i * 1.0 );
+        v2.push_back( i * 1.0 );
+    }
+    BOOST_CHECK( cmp::vector_equal<float>(v1 , v2 ));
+    v1.push_back( 27 );
+    BOOST_CHECK( !cmp::vector_equal<float>(v1 , v2 ));
+    v2.push_back( 27 );
+    BOOST_CHECK( cmp::vector_equal(v1 , v2 ));
+}