From 1ead8dd4589281ac064c3e26b18dca248f012caf Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Mon, 10 Aug 2015 16:05:00 +0200 Subject: [PATCH] Fix bug in floating-point comparison. The function used integer abs() instead of fabs(), in C abs() only takes integers so this is an error. Also needed to add comparison with absolute tolerance to get reasonable behaviour (comparing 0 versus 0.0000000000001 is ok here). --- opm/core/grid/grid.c | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/opm/core/grid/grid.c b/opm/core/grid/grid.c index 26bd8e48..d07017f9 100644 --- a/opm/core/grid/grid.c +++ b/opm/core/grid/grid.c @@ -25,6 +25,7 @@ #include #include #include +#include void @@ -550,14 +551,15 @@ static int memcmp_double(const double * p1 , const double *p2 , size_t num_eleme if (memcmp(p1 , p2 , num_elements * sizeof * p1) == 0) return 0; else { - const double epsilon = 1e-5; + const double abs_epsilon = 1e-8; + const double rel_epsilon = 1e-5; for (size_t index = 0; index < num_elements; index++) { - double diff = abs(p1[index] - p2[index]); - if (diff != 0) { - double sum = abs(p1[index]) + abs(p2[index]); + double diff = fabs(p1[index] - p2[index]); + if (diff > abs_epsilon) { + double sum = fabs(p1[index]) + fabs(p2[index]); - if (diff > sum * epsilon) + if (diff > sum * rel_epsilon) return 1; } }