Diff.cc: merged from develop.

src/Diff.cc

@@ -25,33 +25,38 @@
 #include "pmlist.h"
 #include "cdo_zaxis.h"
 
+static inline void
+diff_kernel(double v1, double v2, size_t &ndiff, bool &dsgn, bool &zero, double &absm, double &relm)
+{
+  const auto absdiff = std::fabs(v1 - v2);
+  if (absdiff > 0.0) ndiff++;
+
+  absm = std::max(absm, absdiff);
+
+  const auto vv = v1 * v2;
+  if (vv < 0.0)
+    dsgn = true;
+  else if (is_equal(vv, 0.0))
+    zero = true;
+  else
+    relm = std::max(relm, absdiff / std::max(std::fabs(v1), std::fabs(v2)));
+}
+
 static void
-diff_kernel(bool hasMissvals, const double v1, const double v2, const double missval1, const double missval2, size_t &ndiff, bool &dsgn, bool &zero,
-            double &absm, double &relm)
+diff_kernel_mv(double v1, double v2, double missval1, double missval2, size_t &ndiff, bool &dsgn, bool &zero, double &absm, double &relm)
 {
   const auto v1isnan = std::isnan(v1);
   const auto v2isnan = std::isnan(v2);
-  const auto v1ismissval = hasMissvals ? DBL_IS_EQUAL(v1, missval1) : false;
-  const auto v2ismissval = hasMissvals ? DBL_IS_EQUAL(v2, missval2) : false;
+  const auto v1ismissval = dbl_is_equal(v1, missval1);
+  const auto v2ismissval = dbl_is_equal(v2, missval2);
   if (v1isnan != v2isnan)
     {
       ndiff++;
       relm = 1.0;
     }
-  else if (!hasMissvals || (!v1ismissval && !v2ismissval))
+  else if (!v1ismissval && !v2ismissval)
     {
-      const double absdiff = std::fabs(v1 - v2);
-      if (absdiff > 0.0) ndiff++;
-
-      absm = std::max(absm, absdiff);
-
-      const auto vv = v1 * v2;
-      if (vv < 0.0)
-        dsgn = true;
-      else if (IS_EQUAL(vv, 0.0))
-        zero = true;
-      else
-        relm = std::max(relm, absdiff / std::max(std::fabs(v1), std::fabs(v2)));
+      diff_kernel(v1, v2, ndiff, dsgn, zero, absm, relm);
     }
   else if (v1ismissval != v2ismissval)
     {
@@ -61,13 +66,17 @@ diff_kernel(bool hasMissvals, const double v1, const double v2, const double mis
 }
 
 static void
-diff(size_t gridsize, bool hasMissvals, const Field_t &field1, const Field_t &field2, size_t &ndiff, bool &dsgn, bool &zero,
-     double &absm, double &relm)
+diff(size_t gridsize, const Field_t &field1, const Field_t &field2, size_t &ndiff, bool &dsgn, bool &zero, double &absm, double &relm)
 {
   auto func = [&](const LambdaField_T_ _field1, const LambdaField_T_ _field2)
     {
-      for (size_t i = 0; i < gridsize; ++i)
-        diff_kernel(hasMissvals, _field1.vec[i], _field2.vec[i], _field1.missval, _field2.missval, ndiff, dsgn, zero, absm, relm);
+      const auto hasMissvals = (_field1.nmiss || _field2.nmiss);
+      if (hasMissvals)
+        for (size_t i = 0; i < gridsize; ++i)
+          diff_kernel_mv(_field1.vec[i], _field2.vec[i], _field1.missval, _field2.missval, ndiff, dsgn, zero, absm, relm);
+      else
+        for (size_t i = 0; i < gridsize; ++i)
+          diff_kernel(_field1.vec[i], _field2.vec[i], ndiff, dsgn, zero, absm, relm);
     };
 
   field_operation2DiffType(func, field1, field2);
@@ -238,12 +247,11 @@ Diff(void *process)
           cdo_read_record(streamID2, field2);
           if (varList2[varID2].nwpv == CDI_COMP) use_real_part(gridsize, field2);
 
-          const auto hasMissvals = (field1.getNumMiss() || field2.getNumMiss());
           size_t ndiff = 0;
           auto dsgn = false, zero = false;
           double absm = 0.0, relm = 0.0;
 
-          diff(gridsize, hasMissvals, field1, field2, ndiff, dsgn, zero, absm, relm);
+          diff(gridsize, field1, field2, ndiff, dsgn, zero, absm, relm);
 
           if (!Options::silentMode || Options::cdoVerbose)
             {