MAINT: add SCIPY_USE_PROPACK env variable (scipy#16361)

* this is effectively a forward port and modernization
of the release branch `PROPACK` shims that were added in
scipygh-15432; in short, `PROPACK` + Windows + some linalg backends
was causing all sorts of trouble, and this has never been resolved

* I've switched to `SCIPY_USE_PROPACK` instead of `USE_PROPACK`
for the opt-in, since this was requested, though the change
between release branches may cause a little confusion (another
release note adjustment to add maybe)

* I think the issues are painful to reproduce; for my part,
I did the following just to check the proper skipping/accounting
of tests:

- `SCIPY_USE_PROPACK=1 python dev.py -j 20 -t scipy/sparse/linalg`
  - `932 passed, 172 skipped, 8 xfailed in 115.57s (0:01:55)`
- `python dev.py -j 20 -t scipy/sparse/linalg`
  - `787 passed, 317 skipped, 8 xfailed in 114.80s (0:01:54)`

* why am I doing this now? well, to be frank the process of manually
backporting this for each release is error-prone, and may cause
additional confusion/debate, which I'd like to avoid. Besides, if it
is broken in `main` we may as well have the shims there as well. I would
point out that you may want to add `SCIPY_USE_PROPACK` to 1 or 2 jobs
in CI? The other reason is that if usage of `PROPACK` spreads, I don't
want to be manually applying more skips/shims on each release (which
I already had to do here with two new tests it seems)

.github/workflows/linux_meson.yml

@@ -101,6 +101,7 @@ jobs:
     - name: Test SciPy
       run: |
         export OMP_NUM_THREADS=2
+        export SCIPY_USE_PROPACK=1
         python dev.py -n -j 2
 
   test_venv_install:

.github/workflows/macos_meson.yml

@@ -112,6 +112,7 @@ jobs:
       run: |
         conda activate scipy-dev
         export OMP_NUM_THREADS=2
+        export SCIPY_USE_PROPACK=1
         python dev.py -n -j 2
 
     - name: Ccache statistics

.github/workflows/windows.yml

@@ -79,6 +79,7 @@ jobs:
       - name: prep-test
         run: |
           echo "PYTHONPATH=${env:installed_path}" >> $env:GITHUB_ENV
+          echo "SCIPY_USE_PROPACK=1" >> $env:GITHUB_ENV
       - name: test
         run: |
           mkdir tmp

scipy/sparse/linalg/_eigen/_svds.py

@@ -1,3 +1,4 @@
+import os
 import numpy as np
 
 from .arpack import _arpack  # type: ignore[attr-defined]
@@ -6,7 +7,11 @@
 from scipy._lib._util import check_random_state
 from scipy.sparse.linalg._interface import LinearOperator, aslinearoperator
 from scipy.sparse.linalg._eigen.lobpcg import lobpcg  # type: ignore[no-redef]
-from scipy.sparse.linalg._svdp import _svdp
+if os.environ.get("SCIPY_USE_PROPACK"):
+    from scipy.sparse.linalg._svdp import _svdp
+    HAS_PROPACK = True
+else:
+    HAS_PROPACK = False
 from scipy.linalg import svd
 
 arpack_int = _arpack.timing.nbx.dtype
@@ -297,6 +302,12 @@ def matmat_XH_X(x):
         eigvec, _ = np.linalg.qr(eigvec)
 
     elif solver == 'propack':
+        if not HAS_PROPACK:
+            raise ValueError("`solver='propack'` is opt-in due "
+                             "to potential issues on Windows, "
+                             "it can be enabled by setting the "
+                             "`SCIPY_USE_PROPACK` environment "
+                             "variable before importing scipy")
         jobu = return_singular_vectors in {True, 'u'}
         jobv = return_singular_vectors in {True, 'vh'}
         irl_mode = (which == 'SM')

scipy/sparse/linalg/_eigen/tests/test_svds.py

@@ -1,3 +1,4 @@
+import os
 import re
 import copy
 import numpy as np
@@ -8,6 +9,10 @@
 from scipy.linalg import hilbert, svd, null_space
 from scipy.sparse import csc_matrix, isspmatrix, spdiags, random
 from scipy.sparse.linalg import LinearOperator, aslinearoperator
+if os.environ.get("SCIPY_USE_PROPACK"):
+    has_propack = True
+else:
+    has_propack = False
 from scipy.sparse.linalg import svds
 from scipy.sparse.linalg._eigen.arpack import ArpackNoConvergence
 
@@ -157,6 +162,8 @@ def test_svds_input_validation_k_1(self, k):
 
         # propack can do complete SVD
         if self.solver == 'propack' and k == 3:
+            if not has_propack:
+                pytest.skip("PROPACK not enabled")
             res = svds(A, k=k, solver=self.solver)
             _check_svds(A, k, *res, check_usvh_A=True, check_svd=True)
             return
@@ -257,6 +264,9 @@ def test_svds_input_validation_return_singular_vectors(self, rsv):
     @pytest.mark.parametrize("k", [3, 5])
     @pytest.mark.parametrize("which", ["LM", "SM"])
     def test_svds_parameter_k_which(self, k, which):
+        if self.solver == 'propack':
+            if not has_propack:
+                pytest.skip("PROPACK not available")
         # check that the `k` parameter sets the number of eigenvalues/
         # eigenvectors returned.
         # Also check that the `which` parameter sets whether the largest or
@@ -274,6 +284,9 @@ def test_svds_parameter_k_which(self, k, which):
 
     # loop instead of parametrize for simplicity
     def test_svds_parameter_tol(self):
+        if self.solver == 'propack':
+            if not has_propack:
+                pytest.skip("PROPACK not available")
         return  # TODO: needs work, disabling for now
         # check the effect of the `tol` parameter on solver accuracy by solving
         # the same problem with varying `tol` and comparing the eigenvalues
@@ -315,6 +328,9 @@ def err(tol):
             assert error > accuracy/10
 
     def test_svd_v0(self):
+        if self.solver == 'propack':
+            if not has_propack:
+                pytest.skip("PROPACK not available")
         # check that the `v0` parameter affects the solution
         n = 100
         k = 1
@@ -347,6 +363,9 @@ def test_svd_v0(self):
             assert_equal(res1a, res1b)
 
     def test_svd_random_state(self):
+        if self.solver == 'propack':
+            if not has_propack:
+                pytest.skip("PROPACK not available")
         # check that the `random_state` parameter affects the solution
         # Admittedly, `n` and `k` are chosen so that all solver pass all
         # these checks. That's a tall order, since LOBPCG doesn't want to
@@ -379,6 +398,10 @@ def test_svd_random_state(self):
                                               np.random.RandomState(0),
                                               np.random.default_rng(0)))
     def test_svd_random_state_2(self, random_state):
+        if self.solver == 'propack':
+            if not has_propack:
+                pytest.skip("PROPACK not available")
+
         n = 100
         k = 1
 
@@ -397,6 +420,10 @@ def test_svd_random_state_2(self, random_state):
                                               np.random.RandomState(0),
                                               np.random.default_rng(0)))
     def test_svd_random_state_3(self, random_state):
+        if self.solver == 'propack':
+            if not has_propack:
+                pytest.skip("PROPACK not available")
+
         n = 100
         k = 5
 
@@ -417,6 +444,9 @@ def test_svd_random_state_3(self, random_state):
     def test_svd_maxiter(self):
         # check that maxiter works as expected: should not return accurate
         # solution after 1 iteration, but should with default `maxiter`
+        if self.solver == 'propack':
+            if not has_propack:
+                pytest.skip("PROPACK not available")
         A = hilbert(6)
         k = 1
         u, s, vh = sorted_svd(A, k)
@@ -443,6 +473,10 @@ def test_svd_maxiter(self):
     @pytest.mark.parametrize("shape", ((5, 7), (6, 6), (7, 5)))
     def test_svd_return_singular_vectors(self, rsv, shape):
         # check that the return_singular_vectors parameter works as expected
+        if self.solver == 'propack':
+            if not has_propack:
+                pytest.skip("PROPACK not available")
+
         rng = np.random.default_rng(0)
         A = rng.random(shape)
         k = 2
@@ -521,6 +555,10 @@ def test_svd_return_singular_vectors(self, rsv, shape):
                                          aslinearoperator))
     def test_svd_simple(self, A, k, real, transpose, lo_type):
 
+        if self.solver == 'propack':
+            if not has_propack:
+                pytest.skip("PROPACK not available")
+
         A = np.asarray(A)
         A = np.real(A) if real else A
         A = A.T if transpose else A
@@ -547,6 +585,9 @@ def test_svd_simple(self, A, k, real, transpose, lo_type):
 
     def test_svd_linop(self):
         solver = self.solver
+        if self.solver == 'propack':
+            if not has_propack:
+                pytest.skip("PROPACK not available")
 
         nmks = [(6, 7, 3),
                 (9, 5, 4),
@@ -720,6 +761,8 @@ def test_svd_LM_zeros_matrix(self, shape, dtype):
     # ARPACK supports only dtype float, complex, or np.float32
     @pytest.mark.parametrize("dtype", (float, complex, np.float32))
     def test_small_sigma(self, shape, dtype):
+        if not has_propack:
+            pytest.skip("PROPACK not enabled")
         # https://github.com/scipy/scipy/pull/11829
         if dtype == complex and self.solver == 'propack':
             pytest.skip("PROPACK unsupported for complex dtype")
@@ -743,6 +786,8 @@ def test_small_sigma(self, shape, dtype):
     @pytest.mark.filterwarnings("ignore:The problem size")
     @pytest.mark.parametrize("dtype", (float, complex, np.float32))
     def test_small_sigma2(self, dtype):
+        if not has_propack:
+            pytest.skip("PROPACK not enabled")
         # https://github.com/scipy/scipy/issues/11406
         if dtype == complex and self.solver == 'propack':
             pytest.skip("PROPACK unsupported for complex dtype")