varrayMinMaxSum(): refactored.

src/Info.cc

@@ -40,9 +40,19 @@ void
 fieldMinMaxSum(const Field &field, double &min, double &max, double &sum)
 {
   if (field.memType == MemType::Float)
-    varrayMinMaxSum(field.size, field.vec_f, min, max, sum);
+    {
+      auto mms = varrayMinMaxSum(field.size, field.vec_f, MinMaxSum(min, max, sum));
+      min = mms.min;
+      max = mms.max;
+      sum = mms.sum;
+    }
   else
-    varrayMinMaxSum(field.size, field.vec_d, min, max, sum);
+    {
+      auto mms = varrayMinMaxSum(field.size, field.vec_d, MinMaxSum(min, max, sum));
+      min = mms.min;
+      max = mms.max;
+      sum = mms.sum;
+    }
 }
 
 size_t

src/array.cc

@@ -80,18 +80,14 @@ varrayMinMaxMV(const size_t len, const Varray<double> &array, const double missv
 }
 
 template <typename T>
-void
-varrayMinMaxSum(const size_t len, const Varray<T> &array, double &rmin, double &rmax, double &rsum)
+MinMaxSum
+varrayMinMaxSum(const size_t len, const Varray<T> &array, MinMaxSum mmsinit)
 {
   // rmin, rmax and rsum will be initialized in Info
+  T rmin = mmsinit.min;
+  T rmax = mmsinit.max;
+  auto rsum = mmsinit.sum;
 
-  /* 
-  gnu 9.2.0: correct result
-  intel 18.0.5: wrong result
-#ifdef HAVE_OPENMP4
-#pragma omp simd reduction(+:rsum) reduction(min:rmin) reduction(max:rmax)
-#endif
-  */
   for (size_t i = 0; i < len; ++i)
     {
       const auto v = array[i];
@@ -99,11 +95,13 @@ varrayMinMaxSum(const size_t len, const Varray<T> &array, double &rmin, double &
       if (v > rmax) rmax = v;
       rsum += v;
     }
+
+  return MinMaxSum(rmin, rmax, rsum, len);
 }
 
 // Explicit instantiation
-template void varrayMinMaxSum(const size_t len, const Varray<float> &array, double &rmin, double &rmax, double &rsum);
-template void varrayMinMaxSum(const size_t len, const Varray<double> &array, double &rmin, double &rmax, double &rsum);
+template MinMaxSum varrayMinMaxSum(const size_t len, const Varray<float> &array, MinMaxSum mmsinit);
+template MinMaxSum varrayMinMaxSum(const size_t len, const Varray<double> &array, MinMaxSum mmsinit);
 
 template <typename T>
 size_t
@@ -137,12 +135,10 @@ template size_t varrayMinMaxSumMV(const size_t len, const Varray<double> &array,
 void
 varrayMinMaxMean(const size_t len, const Varray<double> &array, double &rmin, double &rmax, double &rmean)
 {
-  rmin = DBL_MAX;
-  rmax = -DBL_MAX;
-  double sum = 0.0;
-  varrayMinMaxSum(len, array, rmin, rmax, sum);
-
-  rmean = len ? sum / (double) len : 0;
+  auto mms = varrayMinMaxSum(len, array, MinMaxSum());
+  rmin = mms.min;
+  rmax = mms.max;
+  rmean = len ? mms.sum / (double) len : 0;
 }
 
 size_t

src/array.h

@@ -19,9 +19,21 @@
 
 #include <vector>
 #include <cstddef>
+#include <cfloat>
 #include "compare.h"
 
 
+struct MinMaxSum
+{
+  double min;
+  double max;
+  double sum;
+  size_t n;
+  MinMaxSum() : min(DBL_MAX), max(-DBL_MAX), sum(0), n(0) {};
+  MinMaxSum(double rmin, double rmax, double rsum, size_t rn) : min(rmin), max(rmax), sum(rsum), n(rn) {};
+  MinMaxSum(double rmin, double rmax, double rsum) : min(rmin), max(rmax), sum(rsum) {};
+};
+
 using MinMaxVal = std::pair<double, double>;
 
 template <typename T>
@@ -117,7 +129,7 @@ static inline double FSQRTMN(const double x, const double missval1)
 const char *fpe_errstr(int fpeRaised);
 
 template <typename T>
-void varrayMinMaxSum(size_t len, const Varray<T> &array, double &rmin, double &rmax, double &rsum);
+MinMaxSum varrayMinMaxSum(size_t len, const Varray<T> &array, MinMaxSum mmsinit);
 
 template <typename T>
 size_t varrayMinMaxSumMV(size_t len, const Varray<T> &array, T missval, double &rmin, double &rmax, double &rsum);