OilfieldToolsNet/opm-core
4
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Fix bug in floating-point comparison.

Browse Source 
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).
This commit is contained in:
Atgeirr Flø Rasmussen 2015-08-10 16:05:00 +02:00
parent 2b49064b01
commit 1ead8dd458
1 changed files with 7 additions and 5 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

12
opm/core/grid/grid.c
View File

@ -25,6 +25,7 @@
#include <stdio.h> #include <stdio.h>
#include <stdlib.h> #include <stdlib.h>
#include <string.h> #include <string.h>
#include <math.h>
void 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) if (memcmp(p1 , p2 , num_elements * sizeof * p1) == 0)
return 0; return 0;
else { 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++) { for (size_t index = 0; index < num_elements; index++) {
double diff = abs(p1[index] - p2[index]); double diff = fabs(p1[index] - p2[index]);
if (diff != 0) { if (diff > abs_epsilon) {
double sum = abs(p1[index]) + abs(p2[index]); double sum = fabs(p1[index]) + fabs(p2[index]);
if (diff > sum * epsilon) if (diff > sum * rel_epsilon)
return 1; return 1;
} }
} }