Repository: numpy/numpy
Pull request: ENH: Check for floating point exceptions in dot
Currently dot does not check for floating point exceptions, while similar functions (such as multiply and matmul) do.  Add these checks for consistency.

Note that certain types of floating point exceptions are still not caught:
- Multiplications of 0 with nan or inf are ignored (see #27902)
- Integer overflows are ignored.  These aren't true floating point exceptions, but the behavior is inconsistent with scalar multiply, i.e.
```
>>> np.int16(32000) * np.int16(3)
<python-input-7>:1: RuntimeWarning: overflow encountered in scalar multiply
  np.int16(32000) * np.int16(3)
np.int16(30464)
>>> np.dot(np.int16(32000), np.int16(3))
np.int16(30464)
```

Fixes #14925


Problem statement:
BUG: Dot product does not check for floating point warning flags
The `dot` function (unlike ufuncs, and the ufunc based `matmul`, `@`) does not check for floating point warning/error flags. It would probably better if it would, although it is likely not high priority. The example being:

```
np.dot([1,1], [-np.inf, np.inf])
```
Which should be setting the CPU invalid value produced flag on the NaN result, but it is not caught by NumPy.

(phased out from gh-14788)

Hints from issue discussion:

Commit message: BUG: Check for floating point exceptions in dot

Currently dot does not check for floating point exceptions, while similar
functions (such as multiply and matmul) do.  Add these checks for consistency.

Note that certain types of floating point exceptions are still not caught:
- Multiplications of 0 with nan or inf are ignored (see #27902)
- Integer overflows are ignored.  These aren't true floating point exceptions,
but the behavior is inconsistent with scalar multiply, i.e.
```
>>> np.int16(32000) * np.int16(3)
<python-input-7>:1: RuntimeWarning: overflow encountered in scalar multiply
  np.int16(32000) * np.int16(3)
np.int16(30464)
>>> np.dot(np.int16(32000), np.int16(3))
np.int16(30464)
```

Fixes #14925

Patch to review:
diff --git a/numpy/_core/src/multiarray/multiarraymodule.c b/numpy/_core/src/multiarray/multiarraymodule.c
index d24af97b2fa3..e5b8aeb1ecee 100644
--- a/numpy/_core/src/multiarray/multiarraymodule.c
+++ b/numpy/_core/src/multiarray/multiarraymodule.c
@@ -958,6 +958,20 @@ PyArray_MatrixProduct(PyObject *op1, PyObject *op2)
     return PyArray_MatrixProduct2(op1, op2, NULL);
 }
 
+static int
+check_fperr(const char *func_name)
+{
+    int errmask = 0;
+    if (_get_bufsize_errmask(NULL, &errmask) < 0) {
+        return -1;
+    }
+    if (_check_ufunc_fperr(errmask, func_name) < 0) {
+        return -1;
+    }
+
+    return 0;
+}
+
 /*NUMPY_API
  * Numeric.matrixproduct2(a,v,out)
  * just like inner product but does the swapaxes stuff on the fly
@@ -1008,11 +1022,17 @@ PyArray_MatrixProduct2(PyObject *op1, PyObject *op2, PyArrayObject* out)
         return NULL;
     }
 
+    npy_clear_floatstatus();
+
 #if defined(HAVE_CBLAS)
     if (PyArray_NDIM(ap1) <= 2 && PyArray_NDIM(ap2) <= 2 &&
             (NPY_DOUBLE == typenum || NPY_CDOUBLE == typenum ||
              NPY_FLOAT == typenum || NPY_CFLOAT == typenum)) {
-        return cblas_matrixproduct(typenum, ap1, ap2, out);
+        PyObject *res = cblas_matrixproduct(typenum, ap1, ap2, out);
+        if (check_fperr("dot") < 0) {
+            return NULL;
+        }
+        return res;
     }
 #endif
 
@@ -1102,6 +1122,9 @@ PyArray_MatrixProduct2(PyObject *op1, PyObject *op2, PyArrayObject* out)
         /* only for OBJECT arrays */
         goto fail;
     }
+    if (check_fperr("dot") < 0) {
+        goto fail;
+    }
     Py_DECREF(ap1);
     Py_DECREF(ap2);
 
diff --git a/numpy/_core/tests/test_multiarray.py b/numpy/_core/tests/test_multiarray.py
index 3de42ba8ed80..3cf783c91d95 100644
--- a/numpy/_core/tests/test_multiarray.py
+++ b/numpy/_core/tests/test_multiarray.py
@@ -3345,6 +3345,20 @@ def test_dot(self):
         a.dot(b=b, out=c)
         assert_equal(c, np.dot(a, b))
 
+    @pytest.mark.skipif(IS_WASM, reason="no wasm fp exception support")
+    def test_dot_errstate(self):
+        a = np.array([1, 1])
+        b = np.array([-np.inf, np.inf])
+
+        with np.errstate(invalid='raise'):
+            with pytest.raises(FloatingPointError,
+                    match="invalid value encountered in dot"):
+                np.dot(a, b)
+
+            with pytest.raises(FloatingPointError,
+                    match="invalid value encountered in dot"):
+                np.dot(a[np.newaxis, np.newaxis, ...], b[np.newaxis, ..., np.newaxis])
+
     def test_dot_type_mismatch(self):
         c = 1.
         A = np.array((1, 1), dtype='i,i')


Code review comments:
diff_hunk: @@ -1008,11 +1022,17 @@ PyArray_MatrixProduct2(PyObject *op1, PyObject *op2, PyArrayObject* out)
         return NULL;
     }
 
+    npy_clear_floatstatus();
path: numpy/_core/src/multiarray/multiarraymodule.c
text: @user1:
I would prefer pairing this with the GIL release (i.e. very targeted around the place of interest).

That probably means moving it into `cblas_matrixproduct`.

@author:
I moved this to `cblas_matrixproduct` but I'm not sure if it's quite what you want.  Do you want this immediately before (or after?) the GIL release, or just closer to where the actual work is being done?

diff_hunk: @@ -958,6 +958,20 @@ PyArray_MatrixProduct(PyObject *op1, PyObject *op2)
     return PyArray_MatrixProduct2(op1, op2, NULL);
 }
 
+static int
+check_fperr(const char *func_name)
+{
+    int errmask = 0;
+    if (_get_bufsize_errmask(NULL, &errmask) < 0) {
+        return -1;
+    }
+    if (_check_ufunc_fperr(errmask, func_name) < 0) {
+        return -1;
+    }
+
+    return 0;
+}
path: numpy/_core/src/multiarray/multiarraymodule.c
text: @user1:
Please use `PyUFunc_GiveFloatingpointErrors`, you can find the rest of the pattern elsewhere (e.g. in scalar math).

diff_hunk: @@ -1008,11 +1022,17 @@ PyArray_MatrixProduct2(PyObject *op1, PyObject *op2, PyArrayObject* out)
         return NULL;
     }
 
+    npy_clear_floatstatus();
+
 #if defined(HAVE_CBLAS)
     if (PyArray_NDIM(ap1) <= 2 && PyArray_NDIM(ap2) <= 2 &&
             (NPY_DOUBLE == typenum || NPY_CDOUBLE == typenum ||
              NPY_FLOAT == typenum || NPY_CFLOAT == typenum)) {
-        return cblas_matrixproduct(typenum, ap1, ap2, out);
+        PyObject *res = cblas_matrixproduct(typenum, ap1, ap2, out);
+        if (check_fperr("dot") < 0) {
+            return NULL;
path: numpy/_core/src/multiarray/multiarraymodule.c
text: @user1:
Just to note, this is also incorrect of course:  You would have to throw away `res` here, but you are leaking it.