Merge pull request #5879 from stuartarchibald/fix/5870

Fix erroneous input mutation in linalg routines

numba/np/linalg.py

@@ -782,6 +782,38 @@ def _check_homogeneous_types(func_name, *types):
             raise TypingError(msg, highlighting=False)
 
 
+def _copy_to_fortran_order():
+    pass
+
+
+@overload(_copy_to_fortran_order)
+def ol_copy_to_fortran_order(a):
+    # This function copies the array 'a' into a new array with fortran order.
+    # This exists because the copy routines don't take order flags yet.
+    F_layout = a.layout == 'F'
+    A_layout = a.layout == 'A'
+    def impl(a):
+        if F_layout:
+            # it's F ordered at compile time, just copy
+            acpy = np.copy(a)
+        elif A_layout:
+            # decide based on runtime value
+            flag_f = a.flags.f_contiguous
+            if flag_f:
+                # it's already F ordered, so copy but in a round about way to
+                # ensure that the copy is also F ordered
+                acpy = np.copy(a.T).T
+            else:
+                # it's something else ordered, so let asfortranarray deal with
+                # copying and making it fortran ordered
+                acpy = np.asfortranarray(a)
+        else:
+            # it's C ordered at compile time, asfortranarray it.
+            acpy = np.asfortranarray(a)
+        return acpy
+    return impl
+
+
 @register_jitable
 def _inv_err_handler(r):
     if r != 0:
@@ -810,8 +842,6 @@ def inv_impl(a):
 
     kind = ord(get_blas_kind(a.dtype, "inv"))
 
-    F_layout = a.layout == 'F'
-
     def inv_impl(a):
         n = a.shape[-1]
         if a.shape[-2] != n:
@@ -820,10 +850,7 @@ def inv_impl(a):
 
         _check_finite_matrix(a)
 
-        if F_layout:
-            acpy = np.copy(a)
-        else:
-            acpy = np.asfortranarray(a)
+        acpy = _copy_to_fortran_order(a)
 
         if n == 0:
             return acpy
@@ -929,8 +956,6 @@ def eig_impl(a):
     JOBVL = ord('N')
     JOBVR = ord('V')
 
-    F_layout = a.layout == 'F'
-
     def real_eig_impl(a):
         """
         eig() implementation for real arrays.
@@ -942,10 +967,7 @@ def real_eig_impl(a):
 
         _check_finite_matrix(a)
 
-        if F_layout:
-            acpy = np.copy(a)
-        else:
-            acpy = np.asfortranarray(a)
+        acpy = _copy_to_fortran_order(a)
 
         ldvl = 1
         ldvr = n
@@ -1001,10 +1023,7 @@ def cmplx_eig_impl(a):
 
         _check_finite_matrix(a)
 
-        if F_layout:
-            acpy = np.copy(a)
-        else:
-            acpy = np.asfortranarray(a)
+        acpy = _copy_to_fortran_order(a)
 
         ldvl = 1
         ldvr = n
@@ -1052,8 +1071,6 @@ def eigvals_impl(a):
     JOBVL = ord('N')
     JOBVR = ord('N')
 
-    F_layout = a.layout == 'F'
-
     def real_eigvals_impl(a):
         """
         eigvals() implementation for real arrays.
@@ -1065,10 +1082,7 @@ def real_eigvals_impl(a):
 
         _check_finite_matrix(a)
 
-        if F_layout:
-            acpy = np.copy(a)
-        else:
-            acpy = np.asfortranarray(a)
+        acpy = _copy_to_fortran_order(a)
 
         ldvl = 1
         ldvr = 1
@@ -1127,10 +1141,7 @@ def cmplx_eigvals_impl(a):
 
         _check_finite_matrix(a)
 
-        if F_layout:
-            acpy = np.copy(a)
-        else:
-            acpy = np.asfortranarray(a)
+        acpy = _copy_to_fortran_order(a)
 
         ldvl = 1
         ldvr = 1
@@ -1171,8 +1182,6 @@ def eigh_impl(a):
 
     _check_linalg_matrix(a, "eigh")
 
-    F_layout = a.layout == 'F'
-
     # convert typing floats to numpy floats for use in the impl
     w_type = getattr(a.dtype, "underlying_float", a.dtype)
     w_dtype = np_support.as_dtype(w_type)
@@ -1193,10 +1202,7 @@ def eigh_impl(a):
 
         _check_finite_matrix(a)
 
-        if F_layout:
-            acpy = np.copy(a)
-        else:
-            acpy = np.asfortranarray(a)
+        acpy = _copy_to_fortran_order(a)
 
         w = np.empty(n, dtype=w_dtype)
 
@@ -1225,8 +1231,6 @@ def eigvalsh_impl(a):
 
     _check_linalg_matrix(a, "eigvalsh")
 
-    F_layout = a.layout == 'F'
-
     # convert typing floats to numpy floats for use in the impl
     w_type = getattr(a.dtype, "underlying_float", a.dtype)
     w_dtype = np_support.as_dtype(w_type)
@@ -1247,10 +1251,7 @@ def eigvalsh_impl(a):
 
         _check_finite_matrix(a)
 
-        if F_layout:
-            acpy = np.copy(a)
-        else:
-            acpy = np.asfortranarray(a)
+        acpy = _copy_to_fortran_order(a)
 
         w = np.empty(n, dtype=w_dtype)
 
@@ -1279,8 +1280,6 @@ def svd_impl(a, full_matrices=1):
 
     _check_linalg_matrix(a, "svd")
 
-    F_layout = a.layout == 'F'
-
     # convert typing floats to numpy floats for use in the impl
     s_type = getattr(a.dtype, "underlying_float", a.dtype)
     s_dtype = np_support.as_dtype(s_type)
@@ -1301,10 +1300,7 @@ def svd_impl(a, full_matrices=1):
 
         _check_finite_matrix(a)
 
-        if F_layout:
-            acpy = np.copy(a)
-        else:
-            acpy = np.asfortranarray(a)
+        acpy = _copy_to_fortran_order(a)
 
         ldu = m
         minmn = min(m, n)
@@ -1361,8 +1357,6 @@ def qr_impl(a):
 
     kind = ord(get_blas_kind(a.dtype, "qr"))
 
-    F_layout = a.layout == 'F'
-
     def qr_impl(a):
         n = a.shape[-1]
         m = a.shape[-2]
@@ -1373,10 +1367,7 @@ def qr_impl(a):
         _check_finite_matrix(a)
 
         # copy A as it will be destroyed
-        if F_layout:
-            q = np.copy(a)
-        else:
-            q = np.asfortranarray(a)
+        q = _copy_to_fortran_order(a)
 
         lda = m
 
@@ -1598,9 +1589,6 @@ def lstsq_impl(a, b, rcond=-1.0):
     # B can be 1D or 2D.
     _check_linalg_1_or_2d_matrix(b, "lstsq")
 
-    a_F_layout = a.layout == 'F'
-    b_F_layout = b.layout == 'F'
-
     _check_homogeneous_types("lstsq", a, b)
 
     np_dt = np_support.as_dtype(a.dtype)
@@ -1640,10 +1628,7 @@ def lstsq_impl(a, b, rcond=-1.0):
         maxmn = max(m, n)
 
         # a is destroyed on exit, copy it
-        if a_F_layout:
-            acpy = np.copy(a)
-        else:
-            acpy = np.asfortranarray(a)
+        acpy = _copy_to_fortran_order(a)
 
         # b is overwritten on exit with the solution, copy allocate
         bcpy = np.empty((nrhs, maxmn), dtype=np_dt).T
@@ -1717,9 +1702,6 @@ def solve_impl(a, b):
     _check_linalg_matrix(a, "solve")
     _check_linalg_1_or_2d_matrix(b, "solve")
 
-    a_F_layout = a.layout == 'F'
-    b_F_layout = b.layout == 'F'
-
     _check_homogeneous_types("solve", a, b)
 
     np_dt = np_support.as_dtype(a.dtype)
@@ -1742,10 +1724,7 @@ def solve_impl(a, b):
         _system_check_dimensionally_valid(a, b)
 
         # a is destroyed on exit, copy it
-        if a_F_layout:
-            acpy = np.copy(a)
-        else:
-            acpy = np.asfortranarray(a)
+        acpy = _copy_to_fortran_order(a)
 
         # b is overwritten on exit with the solution, copy allocate
         bcpy = np.empty((nrhs, n), dtype=np_dt).T
@@ -1791,8 +1770,6 @@ def pinv_impl(a, rcond=1.e-15):
 
     numba_xxgemm = _BLAS().numba_xxgemm(a.dtype)
 
-    F_layout = a.layout == 'F'
-
     kind = ord(get_blas_kind(a.dtype, "pinv"))
     JOB = ord('S')
 
@@ -1849,10 +1826,7 @@ def pinv_impl(a, rcond=1.e-15):
 
         _check_finite_matrix(a)
 
-        if F_layout:
-            acpy = np.copy(a)
-        else:
-            acpy = np.asfortranarray(a)
+        acpy = _copy_to_fortran_order(a)
 
         if m == 0 or n == 0:
             return acpy.T.ravel().reshape(a.shape).T
@@ -1992,8 +1966,6 @@ def slogdet_impl(a):
 
     kind = ord(get_blas_kind(a.dtype, "slogdet"))
 
-    F_layout = a.layout == 'F'
-
     diag_walker = _get_slogdet_diag_walker(a)
 
     ONE = a.dtype(1)
@@ -2010,10 +1982,7 @@ def slogdet_impl(a):
 
         _check_finite_matrix(a)
 
-        if F_layout:
-            acpy = np.copy(a)
-        else:
-            acpy = np.asfortranarray(a)
+        acpy = _copy_to_fortran_order(a)
 
         ipiv = np.empty(n, dtype=F_INT_nptype)
 
@@ -2088,8 +2057,6 @@ def _compute_singular_values_impl(a):
     u = np.empty((1, 1), dtype=np_dtype)
     vt = np.empty((1, 1), dtype=np_dtype)
 
-    F_layout = a.layout == 'F'
-
     def sv_function(a):
         """
         Computes singular values.
@@ -2111,10 +2078,7 @@ def sv_function(a):
         ucol = 1
         ldvt = 1
 
-        if F_layout:
-            acpy = np.copy(a)
-        else:
-            acpy = np.asfortranarray(a)
+        acpy = _copy_to_fortran_order(a)
 
         # u and vt are not referenced however need to be
         # allocated (as done above) for MKL as it