ENH: Vectorize float min/max operation with sse2

Improves performance by ~1.5/3.0 for float/double.

3 files changed, 134 insertions, 21 deletions

diff --git a/numpy/core/src/umath/loops.c.src b/numpy/core/src/umath/loops.c.src
index c0287b8c8..db08ad8e2 100644
--- a/numpy/core/src/umath/loops.c.src
+++ b/numpy/core/src/umath/loops.c.src
@@ -1418,6 +1418,9 @@ NPY_NO_EXPORT void
 {
     /*  */
     if (IS_BINARY_REDUCE) {
+        if (run_unary_reduce_simd_@kind@_@TYPE@(args, dimensions, steps)) {
+            return;
+        }
         BINARY_REDUCE_LOOP(@type@) {
             const @type@ in2 = *(@type@ *)ip2;
             io1 = (io1 @OP@ in2 || npy_isnan(io1)) ? io1 : in2;
diff --git a/numpy/core/src/umath/simd.inc.src b/numpy/core/src/umath/simd.inc.src
index 916473a0b..9dad35a52 100644
--- a/numpy/core/src/umath/simd.inc.src
+++ b/numpy/core/src/umath/simd.inc.src
@@ -17,9 +17,14 @@
 
 #include "lowlevel_strided_loops.h"
 #include "npy_config.h"
+/* for NO_FLOATING_POINT_SUPPORT */
+#include "numpy/ufuncobject.h"
 #include <assert.h>
 #include <stdlib.h>
 
+int PyUFunc_getfperr(void);
+void PyUFunc_clearfperr(void);
+
 /*
  * stride is equal to element size and input and destination are equal or
  * don't overlap within one register
@@ -29,21 +34,23 @@
      (npy_is_aligned(args[0], esize) && npy_is_aligned(args[1], esize)) && \
      ((abs(args[1] - args[0]) >= (vsize)) || ((abs(args[1] - args[0]) == 0))))
 
+#define IS_BLOCKABLE_REDUCE(esize, vsize) \
+    (steps[1] == (esize) && abs(args[1] - args[0]) >= (vsize))
+
 
 /* align var to alignment */
-#define UNARY_LOOP_BLOCK_ALIGN_VAR(var, type, alignment)\
+#define LOOP_BLOCK_ALIGN_VAR(var, type, alignment)\
     npy_intp i, peel = npy_aligned_block_offset(var, sizeof(type),\
                                                 alignment, n);\
     for(i = 0; i < peel; i++)
 
-#define UNARY_LOOP_BLOCKED(type, vsize)\
+#define LOOP_BLOCKED(type, vsize)\
     for(; i < npy_blocked_end(peel, sizeof(type), vsize, n);\
             i += (vsize / sizeof(type)))
 
-#define UNARY_LOOP_BLOCKED_END\
+#define LOOP_BLOCKED_END\
     for (; i < n; i++)
 
-
 /*
  * Dispatcher functions
  * decide whether the operation can be vectorized and run it
@@ -58,28 +65,31 @@
  */
 
 /**begin repeat1
- * #func = sqrt, absolute#
+ * #func = sqrt, absolute, minimum, maximum#
+ * #check = IS_BLOCKABLE_UNARY, IS_BLOCKABLE_UNARY, IS_BLOCKABLE_REDUCE, IS_BLOCKABLE_REDUCE#
+ * #name = unary, unary, unary_reduce, unary_reduce#
  */
 
 #if @vector@
 
 /* prototypes */
 static void
-sse2_@func@_@TYPE@(@type@ * op, const @type@ * ip, const npy_intp n);
+sse2_@func@_@TYPE@(@type@ *, @type@ *, const npy_intp n);
 
 #endif
 
 static NPY_INLINE int
-run_unary_simd_@func@_@TYPE@(char **args, npy_intp *dimensions, npy_intp *steps)
+run_@name@_simd_@func@_@TYPE@(char **args, npy_intp *dimensions, npy_intp *steps)
 {
 #if @vector@ && defined HAVE_EMMINTRIN_H
-    if (IS_BLOCKABLE_UNARY(sizeof(@type@), 16)) {
+    if (@check@(sizeof(@type@), 16)) {
         sse2_@func@_@TYPE@((@type@*)args[1], (@type@*)args[0], dimensions[0]);
         return 1;
     }
 #endif
     return 0;
 }
+
 /**end repeat1**/
 
 /**end repeat**/
@@ -89,6 +99,34 @@ run_unary_simd_@func@_@TYPE@(char **args, npy_intp *dimensions, npy_intp *steps)
  * Vectorized operations
  */
 
+#ifdef HAVE_EMMINTRIN_H
+#include <emmintrin.h>
+
+/**begin repeat
+* horizontal reductions on a vector
+* # VOP = min, max#
+*/
+
+static NPY_INLINE npy_float sse2_horizontal_@VOP@___m128(__m128 v)
+{
+    npy_float r;
+    __m128 tmp = _mm_movehl_ps(v, v);                   /* c     d     ... */
+    __m128 m = _mm_@VOP@_ps(v, tmp);                    /* m(ac) m(bd) ... */
+    tmp = _mm_shuffle_ps(m, m, _MM_SHUFFLE(1, 1, 1, 1));/* m(bd) m(bd) ... */
+    _mm_store_ss(&r, _mm_@VOP@_ps(tmp, m));             /* m(acbd) ... */
+    return r;
+}
+
+static NPY_INLINE npy_double sse2_horizontal_@VOP@___m128d(__m128d v)
+{
+    npy_double r;
+    __m128d tmp = _mm_unpackhi_pd(v, v);    /* b     b */
+    _mm_store_sd(&r, _mm_@VOP@_pd(tmp, v)); /* m(ab) m(bb) */
+    return r;
+}
+
+/**end repeat**/
+
 /**begin repeat
  *  #type = npy_float, npy_double#
  *  #TYPE = FLOAT, DOUBLE#
@@ -97,40 +135,38 @@ run_unary_simd_@func@_@TYPE@(char **args, npy_intp *dimensions, npy_intp *steps)
  *  #vtype = __m128, __m128d#
  *  #vpre = _mm, _mm#
  *  #vsuf = ps, pd#
+ *  #nan = NPY_NANF, NPY_NAN#
  */
 
-#ifdef HAVE_EMMINTRIN_H
-#include <emmintrin.h>
-
 
 static void
-sse2_sqrt_@TYPE@(@type@ * op, const @type@ * ip, const npy_intp n)
+sse2_sqrt_@TYPE@(@type@ * op, @type@ * ip, const npy_intp n)
 {
     /* align output to 16 bytes */
-    UNARY_LOOP_BLOCK_ALIGN_VAR(op, @type@, 16) {
+    LOOP_BLOCK_ALIGN_VAR(op, @type@, 16) {
         op[i] = @scalarf@(ip[i]);
     }
     assert(npy_is_aligned(&op[i], 16));
     if (npy_is_aligned(&ip[i], 16)) {
-        UNARY_LOOP_BLOCKED(@type@, 16) {
+        LOOP_BLOCKED(@type@, 16) {
             @vtype@ d = @vpre@_load_@vsuf@(&ip[i]);
             @vpre@_store_@vsuf@(&op[i], @vpre@_sqrt_@vsuf@(d));
         }
     }
     else {
-        UNARY_LOOP_BLOCKED(@type@, 16) {
+        LOOP_BLOCKED(@type@, 16) {
             @vtype@ d = @vpre@_loadu_@vsuf@(&ip[i]);
             @vpre@_store_@vsuf@(&op[i], @vpre@_sqrt_@vsuf@(d));
         }
     }
-    UNARY_LOOP_BLOCKED_END {
+    LOOP_BLOCKED_END {
         op[i] = @scalarf@(ip[i]);
     }
 }
 
 
 static void
-sse2_absolute_@TYPE@(@type@ * op, const @type@ * ip, const npy_intp n)
+sse2_absolute_@TYPE@(@type@ * op, @type@ * ip, const npy_intp n)
 {
     /*
      * get 0x7FFFFFFF mask (everything but signbit set)
@@ -140,34 +176,87 @@ sse2_absolute_@TYPE@(@type@ * op, const @type@ * ip, const npy_intp n)
     const @vtype@ mask = @vpre@_set1_@vsuf@(-0.@c@);
 
     /* align output to 16 bytes */
-    UNARY_LOOP_BLOCK_ALIGN_VAR(op, @type@, 16) {
+    LOOP_BLOCK_ALIGN_VAR(op, @type@, 16) {
         const @type@ tmp = ip[i] > 0 ? ip[i]: -ip[i];
         /* add 0 to clear -0.0 */
         op[i] = tmp + 0;
     }
     assert(npy_is_aligned(&op[i], 16));
     if (npy_is_aligned(&ip[i], 16)) {
-        UNARY_LOOP_BLOCKED(@type@, 16) {
+        LOOP_BLOCKED(@type@, 16) {
             @vtype@ a = @vpre@_load_@vsuf@(&ip[i]);
             @vpre@_store_@vsuf@(&op[i], @vpre@_andnot_@vsuf@(mask, a));
         }
     }
     else {
-        UNARY_LOOP_BLOCKED(@type@, 16) {
+        LOOP_BLOCKED(@type@, 16) {
             @vtype@ a = @vpre@_loadu_@vsuf@(&ip[i]);
             @vpre@_store_@vsuf@(&op[i], @vpre@_andnot_@vsuf@(mask, a));
         }
     }
-    UNARY_LOOP_BLOCKED_END {
+    LOOP_BLOCKED_END {
         const @type@ tmp = ip[i] > 0 ? ip[i]: -ip[i];
         /* add 0 to clear -0.0 */
         op[i] = tmp + 0;
     }
 }
 
+
+/**begin repeat1
+ * #kind = maximum, minimum#
+ * #VOP = max, min#
+ * #OP = >=, <=#
+ **/
+/* arguments swapped as unary reduce has the swapped compared to unary */
+static void
+sse2_@kind@_@TYPE@(@type@ * ip, @type@ * op, const npy_intp n)
+{
+    LOOP_BLOCK_ALIGN_VAR(ip, @type@, 16) {
+        *op = (*op @OP@ ip[i] || npy_isnan(*op)) ? *op : ip[i];
+    }
+    assert(npy_is_aligned(&ip[i], 16));
+    if (i + 2 * 16 / sizeof(@type@) <= n) {
+        /* load the first elements */
+        @vtype@ c = @vpre@_load_@vsuf@((@type@*)&ip[i]);
+#ifdef NO_FLOATING_POINT_SUPPORT
+        @vtype@ cnan = @vpre@_or_@vsuf@(@vpre@_cmpneq_@vsuf@(c, c), cnan);
+#else
+        /* minps/minpd will set invalid flag if nan is encountered */
+        PyUFunc_clearfperr();
 #endif
+        i += 16 / sizeof(@type@);
+
+        LOOP_BLOCKED(@type@, 16) {
+            @vtype@ v = @vpre@_load_@vsuf@((@type@*)&ip[i]);
+            c = @vpre@_@VOP@_@vsuf@(c, v);
+#ifdef NO_FLOATING_POINT_SUPPORT
+            /* check for nan, breaking the loop makes non nan case slow */
+            cnan = @vpre@_or_@vsuf@(@vpre@_cmpneq_@vsuf@(v, v), cnan);
+        }
 
+        if (@vpre@_movemask_@vsuf@(cnan)) {
+            *op = @nan@;
+            return;
+        }
+#else
+        }
+#endif
+        {
+            @type@ tmp = sse2_horizontal_@VOP@_@vtype@(c);
+            if (PyUFunc_getfperr() & UFUNC_FPE_INVALID)
+                *op = @nan@;
+            else
+                *op  = (*op @OP@ tmp || npy_isnan(*op)) ? *op : tmp;
+        }
+    }
+    LOOP_BLOCKED_END {
+        *op  = (*op @OP@ ip[i] || npy_isnan(*op)) ? *op : ip[i];
+    }
+}
+/**end repeat1**/
 
 /**end repeat**/
 
+#endif /* HAVE_EMMINTRIN_H */
+
 #endif
diff --git a/numpy/core/tests/test_umath.py b/numpy/core/tests/test_umath.py
index 6bbb15e6b..c58a0d3f5 100644
--- a/numpy/core/tests/test_umath.py
+++ b/numpy/core/tests/test_umath.py
@@ -687,6 +687,27 @@ class TestSign(TestCase):
             np.seterr(**olderr)
 
 
+class TestMinMax(TestCase):
+    def test_minmax_blocked(self):
+        "simd tests on max/min"
+        for dt in [np.float32, np.float64]:
+            for out, inp, msg in gen_alignment_data(dtype=dt, type='unary',
+                                                    max_size=17):
+                for i in range(inp.size):
+                    inp[:] = np.arange(inp.size, dtype=dt)
+                    inp[i] = np.nan
+                    self.assertTrue(np.isnan(inp.max()),
+                                    msg=repr(inp) + '\n' + msg)
+                    self.assertTrue(np.isnan(inp.min()),
+                                    msg=repr(inp) + '\n' + msg)
+
+                    inp[i] = 1e10
+                    assert_equal(inp.max(), 1e10, err_msg=msg)
+                    inp[i] = -1e10
+                    assert_equal(inp.min(), -1e10, err_msg=msg)
+
+
+
 class TestAbsolute(TestCase):
     def test_abs_blocked(self):
         "simd tests on abs"