diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index f221a492b..eac66f9fe 100644 --- a/CMakeLists_files.cmake +++ b/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 diff --git a/opm/common/data/SimulationDataContainer.cpp b/opm/common/data/SimulationDataContainer.cpp index 9cabf925b..0d7230fd2 100644 --- a/opm/common/data/SimulationDataContainer.cpp +++ b/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( 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( data , other_data )) return false; } else return false; diff --git a/opm/common/util/numeric/cmp.cpp b/opm/common/util/numeric/cmp.cpp deleted file mode 100644 index 7d77aa24b..000000000 --- a/opm/common/util/numeric/cmp.cpp +++ /dev/null @@ -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 . -*/ - -#include -#include - -#include -#include -#include - -#include - - -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& v1, const std::vector& 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& v1, const std::vector& 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 ); - } - - } -} diff --git a/opm/common/util/numeric/cmp.hpp b/opm/common/util/numeric/cmp.hpp index 127f0e4be..95748d239 100644 --- a/opm/common/util/numeric/cmp.hpp +++ b/opm/common/util/numeric/cmp.hpp @@ -22,6 +22,8 @@ #include #include +#include +#include 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() : Compare two values. /// - /// double_ptr_equal(): This compares all the element in the - /// two double * pointers. + /// ptr_equal(): This compares all the element in the + /// two T * pointers. /// - /// double_vector_equal(): This compares all the elements in - /// two std::vector instances. + /// vector_equal(): This compares all the elements in + /// two std::vector instances. /// - /// For both double_vector_equal() and double_ptr_equal() the - /// actual comparison is based on the double_equal() + /// For both vector_equal() and ptr_equal() the + /// actual comparison is based on the scalar_equal() /// 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 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 + bool scalar_equal(T value1, T value2, T abs_eps , T rel_eps) { + static_assert(std::is_floating_point::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 double_equal(double value1, double value2); + bool equal = true; + T diff = fabs(value1 - value2); + if (diff > abs_eps) { + T scale = std::max(fabs(value1), fabs(value2)); - bool double_vector_equal(const std::vector& v1, const std::vector& v2, double abs_eps, double rel_eps); - bool double_vector_equal(const std::vector& v1, const std::vector& v2); + if (diff > scale * rel_eps) { + 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 + bool scalar_equal(T value1, T value2) { + return scalar_equal( value1 , value2 , default_abs_epsilon , default_rel_epsilon ); + } + + template + bool vector_equal(const std::vector& v1, const std::vector& 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( v1[i] , v2[i])) + return false; + } + + return true; + } + + template + bool vector_equal(const std::vector& v1, const std::vector& v2) { + return vector_equal(v1, v2, default_abs_epsilon, default_rel_epsilon); + } + + + template + 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( p1[index] , p2[index] , abs_eps , rel_eps)) { + return false; + } + } + } + return true; + } + + + template + bool ptr_equal(const T* p1, const T* p2, size_t num_elements) { + return ptr_equal(p1, p2, num_elements , default_abs_epsilon, default_rel_epsilon); + } } } diff --git a/tests/test_cmp.cpp b/tests/test_cmp.cpp index e93647437..cc9214e6f 100644 --- a/tests/test_cmp.cpp +++ b/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_EQUAL( false , cmp::double_equal(1,0)); - BOOST_CHECK_EQUAL( false , cmp::double_equal(0,1)); - BOOST_CHECK_EQUAL( false , cmp::double_equal(-1,1)); + BOOST_CHECK( cmp::scalar_equal(1,1)); + BOOST_CHECK_EQUAL( false , cmp::scalar_equal(1,0)); + BOOST_CHECK_EQUAL( false , cmp::scalar_equal(0,1)); + BOOST_CHECK_EQUAL( false , cmp::scalar_equal(-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( v1 , v2)); v1 = 1e-12; v2 = v1 + 0.5*abs_epsilon; - BOOST_CHECK( cmp::double_equal( v1 , v2)); + BOOST_CHECK( cmp::scalar_equal( v1 , v2)); v1 = 7.0; v2 = 7.0; - BOOST_CHECK( cmp::double_equal( v1 , v2)); + BOOST_CHECK( cmp::scalar_equal( v1 , v2)); v1 = -7.0; v2 = -7.0; - BOOST_CHECK( cmp::double_equal( v1 , v2)); + BOOST_CHECK( cmp::scalar_equal( v1 , v2)); v1 = 0; v2 = 0.5 * abs_epsilon; - BOOST_CHECK( cmp::double_equal( v1 , v2)); + BOOST_CHECK( cmp::scalar_equal( v1 , v2)); v1 = 1e7; v2 = 1e7 + 2*abs_epsilon; - BOOST_CHECK( cmp::double_equal( v1 , v2 )); + BOOST_CHECK( cmp::scalar_equal( v1 , v2 )); v1 = 1e7*(1 - abs_epsilon); v2 = 1e7*(1 + rel_epsilon); - BOOST_CHECK( !cmp::double_equal( v1 , v2 )); + BOOST_CHECK( !cmp::scalar_equal( v1 , v2 )); v1 = 1e7*(1 + abs_epsilon); v2 = 1e7*(1 + rel_epsilon); - BOOST_CHECK( cmp::double_equal( v1 , v2 )); + BOOST_CHECK( cmp::scalar_equal( v1 , v2 )); } /* Should be different: */ { v1 = 0; v2 = 1.5 * abs_epsilon; - BOOST_CHECK( !cmp::double_equal( v1 , v2 )); + BOOST_CHECK( !cmp::scalar_equal( v1 , v2 )); v1 = 1e-8; v2 = v1 + 1.5*abs_epsilon; - BOOST_CHECK( !cmp::double_equal( v1 , v2 )); + BOOST_CHECK( !cmp::scalar_equal( v1 , v2 )); v1 = 1; v2 = v1*(1 + 2*rel_epsilon + abs_epsilon); - BOOST_CHECK( !cmp::double_equal( v1 , v2 )); + BOOST_CHECK( !cmp::scalar_equal( v1 , v2 )); v1 = 10; v2 = v1*(1 + 2*rel_epsilon + abs_epsilon); - BOOST_CHECK( !cmp::double_equal( v1 , v2 )); + BOOST_CHECK( !cmp::scalar_equal( v1 , v2 )); v1 = 1e7; v2 = 1e7*(1 + 2*rel_epsilon + abs_epsilon); - BOOST_CHECK( !cmp::double_equal( v1 , v2 )); + BOOST_CHECK( !cmp::scalar_equal( v1 , v2 )); } } +/* Ensure that float instantiation works. */ +BOOST_AUTO_TEST_CASE(TestFloatcmp) { + std::vector v1; + std::vector 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(v1 , v2 )); + v1.push_back( 27 ); + BOOST_CHECK( !cmp::vector_equal(v1 , v2 )); + v2.push_back( 27 ); + BOOST_CHECK( cmp::vector_equal(v1 , v2 )); +}