Add a generic packet ops corresponding to {std}::fmin and {std}::fmax. The non-sensical NaN-propagation rules for std::min std::max implemented by pmin and pmax in Eigen is a longstanding source og confusion and bug report. This change is a first step towards addressing it, as discussing in issue #564.

Eigen/src/Core/GenericPacketMath.h

@@ -215,11 +215,13 @@ pmul(const bool& a, const bool& b) { return a && b; }
 template<typename Packet> EIGEN_DEVICE_FUNC inline Packet
 pdiv(const Packet& a, const Packet& b) { return a/b; }
 
-/** \internal \returns the min of \a a and \a b  (coeff-wise) */
+/** \internal \returns the min of \a a and \a b  (coeff-wise).
+Equivalent to std::min(a, b), so if either a or b is NaN, a is returned. */
 template<typename Packet> EIGEN_DEVICE_FUNC inline Packet
 pmin(const Packet& a, const Packet& b) { return numext::mini(a, b); }
 
-/** \internal \returns the max of \a a and \a b  (coeff-wise) */
+/** \internal \returns the max of \a a and \a b  (coeff-wise)
+Equivalent to std::max(a, b), so if either a or b is NaN, a is returned.*/
 template<typename Packet> EIGEN_DEVICE_FUNC inline Packet
 pmax(const Packet& a, const Packet& b) { return numext::maxi(a, b); }
 
@@ -633,6 +635,23 @@ Packet print(const Packet& a) { using numext::rint; return rint(a); }
 template<typename Packet> EIGEN_DECLARE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS
 Packet pceil(const Packet& a) { using numext::ceil; return ceil(a); }
 
+/** \internal \returns the min of \a a and \a b  (coeff-wise)
+    Equivalent to std::fmin(a, b). Only if both a and b are NaN is NaN returned.
+*/
+template<typename Packet> EIGEN_DEVICE_FUNC inline Packet
+pfmin(const Packet& a, const Packet& b) {
+  Packet not_nan_mask = pcmp_eq(a, a);
+  return pselect(not_nan_mask, pmin(a, b), b);
+}
+
+/** \internal \returns the max of \a a and \a b  (coeff-wise)
+    Equivalent to std::fmax(a, b). Only if both a and b are NaN is NaN returned.*/
+template<typename Packet> EIGEN_DEVICE_FUNC inline Packet
+pfmax(const Packet& a, const Packet& b) {
+  Packet not_nan_mask = pcmp_eq(a, a);
+  return pselect(not_nan_mask, pmax(a, b), b);
+}
+
 /***************************************************************************
 * The following functions might not have to be overwritten for vectorized types
 ***************************************************************************/

test/packetmath.cpp

@@ -783,6 +783,15 @@ void packetmath_notcomplex() {
 
   CHECK_CWISE2_IF(PacketTraits::HasMin, (std::min), internal::pmin);
   CHECK_CWISE2_IF(PacketTraits::HasMax, (std::max), internal::pmax);
+#if EIGEN_HAS_CXX11_MATH
+  using std::fmin;
+  using std::fmax;
+#else
+  using ::fmin;
+  using ::fmax;
+#endif
+  CHECK_CWISE2_IF(PacketTraits::HasMin, fmin, internal::pfmin);
+  CHECK_CWISE2_IF(PacketTraits::HasMax, fmax, internal::pfmax);
   CHECK_CWISE1(numext::abs, internal::pabs);
   CHECK_CWISE2_IF(PacketTraits::HasAbsDiff, REF_ABS_DIFF, internal::pabsdiff);
 
@@ -815,6 +824,17 @@ void packetmath_notcomplex() {
       for (unsigned int i = 0; i < sizeof(Scalar); ++i) data1_bits[k * sizeof(Scalar) + i] = 0x00;
     }
   }
+
+  for (int i = 0; i < PacketSize; ++i) {
+    data1[i] = internal::random<bool>() ? std::numeric_limits<Scalar>::quiet_NaN() : Scalar(0);
+    data1[i + PacketSize] = internal::random<bool>() ? std::numeric_limits<Scalar>::quiet_NaN() : Scalar(0);
+  }
+  // Test NaN propagation for pmin and pmax. It should be equivalent to std::min.
+  CHECK_CWISE2_IF(PacketTraits::HasMin, (std::min), internal::pmin);
+  CHECK_CWISE2_IF(PacketTraits::HasMax, (std::max), internal::pmax);
+  // Test NaN propagation for pfmin and pfmax. It should be equivalent to std::fmin.
+  CHECK_CWISE2_IF(PacketTraits::HasMin, fmin, internal::pfmin);
+  CHECK_CWISE2_IF(PacketTraits::HasMax, fmax, internal::pfmax);
 }
 
 template <>