BUG: core/umath_linalg: ensure FP error status reflects LAPACK error status

1 files changed, 69 insertions, 0 deletions

diff --git a/numpy/core/src/umath/umath_linalg.c.src b/numpy/core/src/umath/umath_linalg.c.src
index a63bcce9a..c9d55f635 100644
--- a/numpy/core/src/umath/umath_linalg.c.src
+++ b/numpy/core/src/umath/umath_linalg.c.src
@@ -374,6 +374,12 @@ offset_ptr(void* ptr, ptrdiff_t offset)
     return (void*)((npy_uint8*)ptr + offset);
 }
 
+static inline void
+clear_fp_errors()
+{
+    PyUFunc_getfperr();
+}
+
 /*
  *****************************************************************************
  **                      Some handy constants                               **
@@ -1780,6 +1786,7 @@ static inline void
     size_t outer_dim = *dimensions++;
     size_t op_count = (JOBZ=='N')?2:3;
     EIGH_PARAMS_t eigh_params;
+    int error_occurred = 0;
 
     for (iter=0; iter < op_count; ++iter) {
         outer_steps[iter] = (ptrdiff_t) steps[iter];
@@ -1824,6 +1831,7 @@ static inline void
                 }
             } else {
                 /* lapack fail, set result to nan */
+                error_occurred = 1;
                 nan_@BASETYPE@_matrix(args[1], &eigenvalues_out_ld);
                 if ('V' == eigh_params.JOBZ) {
                     nan_@TYPE@_matrix(args[2], &eigenvectors_out_ld);
@@ -1834,6 +1842,11 @@ static inline void
 
         release_@lapack_func@(&eigh_params);
     }
+
+    if (!error_occurred) {
+        /* lapack success, clear FP error status */
+        clear_fp_errors();
+    }
 }
 /**end repeat**/
 
@@ -1960,6 +1973,7 @@ static void
 {
     GESV_PARAMS_t params;
     fortran_int n, nrhs;
+    int error_occurred = 0;
     INIT_OUTER_LOOP_3
 
     n = (fortran_int)dimensions[0];
@@ -1979,12 +1993,18 @@ static void
             if (!not_ok) {
                 delinearize_@TYPE@_matrix(args[2], params.B, &r_out);
             } else {
+                error_occurred = 1;
                 nan_@TYPE@_matrix(args[2], &r_out);
             }
         END_OUTER_LOOP
 
         release_@lapack_func@(&params);
     }
+
+    if (!error_occurred) {
+        /* lapack success, clear FP error status */
+        clear_fp_errors();
+    }
 }
 
 static void
@@ -1992,6 +2012,7 @@ static void
               void *NPY_UNUSED(func))
 {
     GESV_PARAMS_t params;
+    int error_occurred;
     fortran_int n;
     INIT_OUTER_LOOP_3
 
@@ -2010,12 +2031,18 @@ static void
             if (!not_ok) {
                 delinearize_@TYPE@_matrix(args[2], params.B, &r_out);
             } else {
+                error_occurred = 1;
                 nan_@TYPE@_matrix(args[2], &r_out);
             }
         END_OUTER_LOOP
 
         release_@lapack_func@(&params);
     }
+
+    if (!error_occurred) {
+        /* lapack success, clear FP error status */
+        clear_fp_errors();
+    }
 }
 
 static void
@@ -2024,6 +2051,7 @@ static void
 {
     GESV_PARAMS_t params;
     fortran_int n;
+    int error_occurred = 0;
     INIT_OUTER_LOOP_2
 
     n = (fortran_int)dimensions[0];
@@ -2040,12 +2068,18 @@ static void
             if (!not_ok) {
                 delinearize_@TYPE@_matrix(args[1], params.B, &r_out);
             } else {
+                error_occurred = 1;
                 nan_@TYPE@_matrix(args[1], &r_out);
             }
         END_OUTER_LOOP
 
         release_@lapack_func@(&params);
     }
+
+    if (!error_occurred) {
+        /* lapack success, clear FP error status */
+        clear_fp_errors();
+    }
 }
 /**end repeat**/
 
@@ -2116,6 +2150,7 @@ static void
 @TYPE@_cholesky(char uplo, char **args, npy_intp *dimensions, npy_intp *steps)
 {
     POTR_PARAMS_t params;
+    int error_occurred = 0;
     fortran_int n;
     INIT_OUTER_LOOP_2
 
@@ -2133,11 +2168,17 @@ static void
                 triu_@TYPE@_matrix(params.A, params.N);
                 delinearize_@TYPE@_matrix(args[1], params.A, &r_out);
             } else {
+                error_occurred = 1;
                 nan_@TYPE@_matrix(args[1], &r_out);
             }
         END_OUTER_LOOP
         release_@lapack_func@(&params);
     }
+
+    if (!error_occurred) {
+        /* lapack success, clear FP error status */
+        clear_fp_errors();
+    }
 }
 
 static void
@@ -2551,6 +2592,7 @@ static inline void
     size_t iter;
     size_t outer_dim = *dimensions++;
     size_t op_count = 2;
+    int error_occurred = 0;
     GEEV_PARAMS_t geev_params;
 
     STACK_TRACE;
@@ -2613,6 +2655,7 @@ static inline void
                                                      &vr_out);
             } else {
                 /* geev failed */
+                error_occurred = 1;
                 nan_@COMPLEXTYPE@_matrix(*arg_iter++, &w_out);
                 if ('V' == geev_params.JOBVL)
                     nan_@COMPLEXTYPE@_matrix(*arg_iter++, &vl_out);
@@ -2624,6 +2667,11 @@ static inline void
 
         release_@lapack_func@(&geev_params);
     }
+
+    if (!error_occurred) {
+        /* lapack success, clear FP error status */
+        clear_fp_errors();
+    }
 }
 
 static void
@@ -2990,6 +3038,7 @@ static inline void
                    npy_intp* steps)
 {
     ptrdiff_t outer_steps[3];
+    int error_occurred = 0;
     size_t iter;
     size_t outer_dim = *dimensions++;
     size_t op_count = (JOBZ=='N')?2:4;
@@ -3046,6 +3095,7 @@ static inline void
                     delinearize_@TYPE@_matrix(args[3], params.VT, &v_out);
                 }
             } else {
+                error_occurred = 1;
                 if ('N' == params.JOBZ) {
                     nan_@REALTYPE@_matrix(args[1], &s_out);
                 } else {
@@ -3059,6 +3109,11 @@ static inline void
 
         release_@lapack_func@(&params);
     }
+
+    if (!error_occurred) {
+        /* lapack success, clear FP error status */
+        clear_fp_errors();
+    }
 }
 /**end repeat*/
 
@@ -3464,6 +3519,7 @@ static void
                 )
 {
     POTRS_PARAMS_t params;
+    int error_occurred = 0;
     fortran_int n, nrhs;
     INIT_OUTER_LOOP_3
 
@@ -3489,12 +3545,18 @@ static void
             if (!not_ok) {
                 delinearize_@TYPE@_matrix(args[2], params.B, &r_out);
             } else {
+                error_occurred = 1;
                 nan_@TYPE@_matrix(args[2], &r_out);
             }
         END_OUTER_LOOP
 
         release_@lapack_func@(&params);
     }
+
+    if (!error_occurred) {
+        /* lapack success, clear FP error status */
+        clear_fp_errors();
+    }
 }
 
 static void
@@ -3617,6 +3679,7 @@ static inline void
              npy_intp *steps)
 {
     POTRI_PARAMS_t params;
+    int error_occurred = 0;
     fortran_int n;
     INIT_OUTER_LOOP_2
 
@@ -3635,12 +3698,18 @@ static inline void
                 make_symmetric_@TYPE@_matrix(params.UPLO, params.N, params.A);
                 delinearize_@TYPE@_matrix(args[1], params.A, &r_out);
             } else {
+                error_occurred = 1;
                 nan_@TYPE@_matrix(args[1], &r_out);
             }
         END_OUTER_LOOP
 
         release_@potri@(&params);
     }
+
+    if (!error_occurred) {
+        /* lapack success, clear FP error status */
+        clear_fp_errors();
+    }
 }
 
 static void