File intc.patch of Package python-scipy

Overview Repositories Revisions Requests Users Attributes Meta

File intc.patch of Package python-scipy (Revision 22428b1546486eae4edaaa0ad313a3cc)

Currently displaying revision 22428b1546486eae4edaaa0ad313a3cc , Show latest

From 8501b7c2fb8a7121aeef94489ece988043c463d0 Mon Sep 17 00:00:00 2001
From: CJ Carey <perimosocordiae@gmail.com>
Date: Sun, 2 Jul 2023 16:29:22 -0400
Subject: [PATCH] BUG: sparse.linalg: Cast index arrays to intc before calling
 SuperLU functions (#18644)

* BUG: sparse.linalg: Cast to intc before SuperLU

---------

Co-authored-by: Stefan van der Walt <stefanv@berkeley.edu>
Co-authored-by: Dan Schult <dschult@colgate.edu>
---
 scipy/sparse/linalg/_dsolve/linsolve.py       | 27 +++++++++++--
 .../linalg/_dsolve/tests/test_linsolve.py     | 38 ++++++++++---------
 2 files changed, 44 insertions(+), 21 deletions(-)

diff --git a/scipy/sparse/linalg/_dsolve/linsolve.py b/scipy/sparse/linalg/_dsolve/linsolve.py
index 0c408f737266..78ba447bcd44 100644
--- a/scipy/sparse/linalg/_dsolve/linsolve.py
+++ b/scipy/sparse/linalg/_dsolve/linsolve.py
@@ -123,6 +123,21 @@ def _get_umf_family(A):
 
     return family, A_new
 
+def _safe_downcast_indices(A):
+    # check for safe downcasting
+    max_value = np.iinfo(np.intc).max
+
+    if A.indptr[-1] > max_value:  # indptr[-1] is max b/c indptr always sorted
+        raise ValueError("indptr values too large for SuperLU")
+
+    if max(*A.shape) > max_value:  # only check large enough arrays
+        if np.any(A.indices > max_value):
+            raise ValueError("indices values too large for SuperLU")
+
+    indices = A.indices.astype(np.intc, copy=False)
+    indptr = A.indptr.astype(np.intc, copy=False)
+    return indices, indptr
+
 def spsolve(A, b, permc_spec=None, use_umfpack=True):
     """Solve the sparse linear system Ax=b, where b may be a vector or a matrix.
 
@@ -269,8 +284,10 @@ def spsolve(A, b, permc_spec=None, use_umfpack=True):
             else:
                 flag = 0  # CSR format
 
+            indices = A.indices.astype(np.intc, copy=False)
+            indptr = A.indptr.astype(np.intc, copy=False)
             options = dict(ColPerm=permc_spec)
-            x, info = _superlu.gssv(N, A.nnz, A.data, A.indices, A.indptr,
+            x, info = _superlu.gssv(N, A.nnz, A.data, indices, indptr,
                                     b, flag, options=options)
             if info != 0:
                 warn("Matrix is exactly singular", MatrixRankWarning)
@@ -402,6 +419,8 @@ def csc_construct_func(*a, cls=type(A)):
     if (M != N):
         raise ValueError("can only factor square matrices")  # is this true?
 
+    indices, indptr = _safe_downcast_indices(A)
+
     _options = dict(DiagPivotThresh=diag_pivot_thresh, ColPerm=permc_spec,
                     PanelSize=panel_size, Relax=relax)
     if options is not None:
@@ -411,7 +430,7 @@ def csc_construct_func(*a, cls=type(A)):
     if (_options["ColPerm"] == "NATURAL"):
         _options["SymmetricMode"] = True
 
-    return _superlu.gstrf(N, A.nnz, A.data, A.indices, A.indptr,
+    return _superlu.gstrf(N, A.nnz, A.data, indices, indptr,
                           csc_construct_func=csc_construct_func,
                           ilu=False, options=_options)
 
@@ -495,6 +514,8 @@ def csc_construct_func(*a, cls=type(A)):
     if (M != N):
         raise ValueError("can only factor square matrices")  # is this true?
 
+    indices, indptr = _safe_downcast_indices(A)
+
     _options = dict(ILU_DropRule=drop_rule, ILU_DropTol=drop_tol,
                     ILU_FillFactor=fill_factor,
                     DiagPivotThresh=diag_pivot_thresh, ColPerm=permc_spec,
@@ -506,7 +527,7 @@ def csc_construct_func(*a, cls=type(A)):
     if (_options["ColPerm"] == "NATURAL"):
         _options["SymmetricMode"] = True
 
-    return _superlu.gstrf(N, A.nnz, A.data, A.indices, A.indptr,
+    return _superlu.gstrf(N, A.nnz, A.data, indices, indptr,
                           csc_construct_func=csc_construct_func,
                           ilu=True, options=_options)
 
diff --git a/scipy/sparse/linalg/_dsolve/tests/test_linsolve.py b/scipy/sparse/linalg/_dsolve/tests/test_linsolve.py
index 4740685a3dd7..00e0d1886dcc 100644
--- a/scipy/sparse/linalg/_dsolve/tests/test_linsolve.py
+++ b/scipy/sparse/linalg/_dsolve/tests/test_linsolve.py
@@ -220,8 +220,11 @@ def test_singular_gh_3312(self):
         except RuntimeError:
             pass
 
-    def test_twodiags(self):
-        A = spdiags([[1, 2, 3, 4, 5], [6, 5, 8, 9, 10]], [0, 1], 5, 5)
+    @pytest.mark.parametrize('format', ['csc', 'csr'])
+    @pytest.mark.parametrize('idx_dtype', [np.int32, np.int64])
+    def test_twodiags(self, format: str, idx_dtype: np.dtype):
+        A = spdiags([[1, 2, 3, 4, 5], [6, 5, 8, 9, 10]], [0, 1], 5, 5,
+                    format=format)
         b = array([1, 2, 3, 4, 5])
 
         # condition number of A
@@ -230,13 +233,12 @@ def test_twodiags(self):
         for t in ['f','d','F','D']:
             eps = finfo(t).eps  # floating point epsilon
             b = b.astype(t)
+            Asp = A.astype(t)
+            Asp.indices = Asp.indices.astype(idx_dtype, copy=False)
+            Asp.indptr = Asp.indptr.astype(idx_dtype, copy=False)
 
-            for format in ['csc','csr']:
-                Asp = A.astype(t).asformat(format)
-
-                x = spsolve(Asp,b)
-
-                assert_(norm(b - Asp@x) < 10 * cond_A * eps)
+            x = spsolve(Asp, b)
+            assert_(norm(b - Asp@x) < 10 * cond_A * eps)
 
     def test_bvector_smoketest(self):
         Adense = array([[0., 1., 1.],
@@ -442,16 +444,18 @@ def setup_method(self):
         n = 40
         d = arange(n) + 1
         self.n = n
-        self.A = spdiags((d, 2*d, d[::-1]), (-3, 0, 5), n, n)
+        self.A = spdiags((d, 2*d, d[::-1]), (-3, 0, 5), n, n, format='csc')
         random.seed(1234)
 
-    def _smoketest(self, spxlu, check, dtype):
+    def _smoketest(self, spxlu, check, dtype, idx_dtype):
         if np.issubdtype(dtype, np.complexfloating):
             A = self.A + 1j*self.A.T
         else:
             A = self.A
 
         A = A.astype(dtype)
+        A.indices = A.indices.astype(idx_dtype, copy=False)
+        A.indptr = A.indptr.astype(idx_dtype, copy=False)
         lu = spxlu(A)
 
         rng = random.RandomState(1234)
@@ -489,10 +493,9 @@ def check(A, b, x, msg=""):
             r = A @ x
             assert_(abs(r - b).max() < 1e3*eps, msg)
 
-        self._smoketest(splu, check, np.float32)
-        self._smoketest(splu, check, np.float64)
-        self._smoketest(splu, check, np.complex64)
-        self._smoketest(splu, check, np.complex128)
+        for dtype in [np.float32, np.float64, np.complex64, np.complex128]:
+            for idx_dtype in [np.int32, np.int64]:
+                self._smoketest(splu, check, dtype, idx_dtype)
 
     @sup_sparse_efficiency
     def test_spilu_smoketest(self):
@@ -508,10 +511,9 @@ def check(A, b, x, msg=""):
             if b.dtype in (np.float64, np.complex128):
                 errors.append(err)
 
-        self._smoketest(spilu, check, np.float32)
-        self._smoketest(spilu, check, np.float64)
-        self._smoketest(spilu, check, np.complex64)
-        self._smoketest(spilu, check, np.complex128)
+        for dtype in [np.float32, np.float64, np.complex64, np.complex128]:
+            for idx_dtype in [np.int32, np.int64]:
+                self._smoketest(spilu, check, dtype, idx_dtype)
 
         assert_(max(errors) > 1e-5)

Places

File intc.patch of Package python-scipy (Revision 22428b1546486eae4edaaa0ad313a3cc)

Places