Feature sparse linalg solvers by abagusetty · Pull Request #2841 · IntelPython/dpnp

abagusetty · 2026-04-09T03:38:05Z

Adds support for from dpnp.scipy.sparse.linalg import LinearOperator, cg, gmres, minres
Fixes: #2831

Have you provided a meaningful PR description?
Have you added a test, reproducer or referred to an issue with a reproducer?
Have you tested your changes locally for CPU and GPU devices?
Have you made sure that new changes do not introduce compiler warnings?
Have you checked performance impact of proposed changes?
Have you added documentation for your changes, if necessary?
Have you added your changes to the changelog?

…oneMKL hooks - _interface.py: add full operator algebra (.H, .T, +, *, **, neg), _AdjointLinearOperator, _TransposedLinearOperator, _SumLinearOperator, _ProductLinearOperator, _ScaledLinearOperator, _PowerLinearOperator, IdentityOperator, MatrixLinearOperator, _AdjointMatrixOperator, _CustomLinearOperator factory dispatch; extend aslinearoperator to handle dpnp sparse and dense arrays - _iterative.py: add _make_system (dtype validation, preconditioner wiring, working dtype selection); add _make_fast_matvec CSR/oneMKL SpMV hook; fix GMRES Arnoldi inner product to single oneMKL BLAS gemv (dpnp.dot) instead of slow Python vdot loop; offload Hessenberg lstsq to numpy.linalg.lstsq (CPU, matches CuPy); fix SciPy host-fallback tol->rtol deprecation via _scipy_tol_kwarg; add preconditioner support to CG; keep MINRES as SciPy-backed stub Refs: CuPy v14.0.1 cupyx/scipy/sparse/linalg/_interface.py, cupyx/scipy/sparse/linalg/_iterative.py"

…gmres, minres Modeled after CuPy's cupyx_tests/scipy_tests/sparse_tests/test_linalg.py. Covers: - LinearOperator: shape, dtype inference, matvec/rmatvec/matmat, subclassing, __matmul__, __call__, edge cases - aslinearoperator: dense array, duck-type, identity passthrough, rmatvec from dense, invalid inputs - cg: SPD convergence, scipy reference match, x0 warm start, b_ndim=2, callback, atol, LinearOperator path, invalid inputs, non-convergence info check - gmres: diag-dominant convergence, scipy reference match, restart variants, x0, b_ndim=2, callbacks, complex systems, atol, non-convergence info check, Hilbert-matrix stress test - minres: SPD, symmetric-indefinite, scipy reference, shift parameter, non-square guard, LinearOperator path, callback - Integration: parametric (n, dtype) cross-solver tests via LinearOperator - Import smoke tests: __all__ completeness

- Use dpnp.tests.helper: assert_dtype_allclose, generate_random_numpy_array, get_all_dtypes, get_float_complex_dtypes, has_support_aspect64 - Use dpnp.tests.third_party.cupy testing harness (with_requires, etc.) - Use numpy.testing assert_allclose / assert_array_equal / assert_raises - Use dpnp.asnumpy() instead of numpy.asarray() - Use pytest parametrize ids matching existing test conventions - Use is_scipy_available() helper from tests/helper.py - Strict class-per-solver organisation matching TestCholesky / TestDet etc.

…or dtype Two bugs fixed: 1. _init_dtype() was calling dpnp.zeros(n) which defaults to float64, so a float32 matvec would upcast and return float64, making the inferred dtype wrong. Fix: use dpnp.zeros(n, dtype=dpnp.int8) as SciPy/CuPy do — any numeric matvec will promote int8 to its own dtype. 2. _CustomLinearOperator.__init__ called _init_dtype() even when an explicit dtype was already supplied, overwriting the caller's value. Fix: _init_dtype() now short-circuits when self.dtype is already set.

…iterative.py

…ption handling Align gemv.cpp with the conventions established in blas/gemm.cpp: Headers added: - ext/common.hpp (dpctl_td_ns, consistent with other extensions) - utils/memory_overlap.hpp (MemoryOverlap guard on x vs y) - utils/output_validation.hpp (CheckWritable + AmpleMemory on y) - utils/type_utils.hpp (validate_type_for_device<T> in impl) - <sstream> (needed for stringstream error_msg) Exception handling added in sparse_gemv_impl(): - try/catch(oneapi::mkl::exception) around all oneMKL sparse calls - try/catch(sycl::exception) around all oneMKL sparse calls - release_matrix_handle cleanup in the exception error path - throw std::runtime_error with descriptive message on catch Input validation added in sparse_gemv(): - ndim checks: x and y must be 1-D - queues_are_compatible() across all 5 USM arrays - MemoryOverlap()(x, y) aliasing guard - CheckWritable::throw_if_not_writable(y) - AmpleMemory::throw_if_not_ample(y, num_rows) - keep_args_alive() at function exit (was missing, returning empty event)

… table Modeled after blas/gemm.cpp (2-D table: value type x index type) and blas/gemv.cpp (dispatch vector pattern with ContigFactory + init_dispatch_table). Changes: - Add sparse/types_matrix.hpp with SparseGemvTypePairSupportFactory<Tv, Ti> encoding the 4 supported combinations: {float32,float64} x {int32,int64} - Rewrite sparse_gemv_impl() to take typeless char* pointers (matching the blas gemv_impl signature style) — type info flows through template params only, no runtime branching inside the impl - Replace the 60-line if/else val_typenum/idx_typenum chain in sparse_gemv() with a 2-D dispatch table lookup (gemv_dispatch_table[val_id][idx_id]) - Rename init_sparse_gemv_dispatch_vector -> init_sparse_gemv_dispatch_table and implement it via init_dispatch_table<> from ext/common.hpp - All validation guards and exception handling from prior commit are preserved

…se_gemv_dispatch_table Follows the rename made in gemv.cpp when the dispatch mechanism was changed from a 1-D vector to a 2-D table (value type x index type). All other declarations (sparse_gemv signature, parameters) are unchanged.

The oneMKL 2025-2 sparse BLAS API deprecated the old 8-argument set_csr_data(queue, handle, nrows, ncols, index_base, row_ptr, col_ind, values, deps) overload in favour of a new signature that takes the sparse matrix handle as `spmat` and adds an explicit `nnz` argument: set_csr_data(queue, spmat, nrows, ncols, nnz, index_base, row_ptr, col_ind, values, deps) Fixes: - Replace old set_csr_data call with the new nnz-aware signature - Silences the resulting -Wunused-parameter warning on `nnz` (now used) - No functional change; all other logic is unchanged

…tring Line 477: `hasattr(A, "rmatmat\")` had a Markdown-escaped backslash leaked into the Python source, causing an unterminated string literal. Fixed to `hasattr(A, "rmatmat")`.

dpnp.ndarray blocks implicit NumPy conversion via __array__ to prevent silent dtype=object arrays. All test assertions must use .asnumpy() to materialize device arrays onto the host explicitly. Also replaces numpy.asarray(x_dp) in _rel_residual helper.

…dation order - _iterative.py: raise NotImplementedError for M != None *before* the _HOST_N_THRESHOLD SciPy fast-path in cg() and gmres(), so the contract is enforced regardless of system size (fixes test_cg_preconditioner_unsupported_raises, test_gmres_preconditioner_unsupported_raises). - _iterative.py: validate callback_type and raise NotImplementedError for 'pr_norm' *before* the _HOST_N_THRESHOLD branch in gmres(), so small-n systems also see the error (fixes test_gmres_callback_type_pr_norm_raises). - _iterative.py: pass callback_type='legacy' to scipy.sparse.linalg.gmres when delegating on the fast path to suppress SciPy DeprecationWarning. - test_scipy_sparse_linalg.py: add dtype=numpy.float64 to expected arange() calls in test_identity_operator and test_gmres_happy_breakdown so strict NumPy 2.0 dtype-equality checks pass (float64 result vs int64 expected).

… port SciPy corner cases

…ecurrence

- Replace .asnumpy() method calls with dpnp.asnumpy() module fn (asnumpy is not an ndarray method in dpnp; it is a top-level fn) - Fix dpnp.any(x) ambiguous truth value in x0 zero-check; replace with explicit `x0 is not None` guard for r0 initialisation - Fix V_mat.T.conj() -> dpnp.conj(V_mat.T) in GMRES Arnoldi step - Guard minres beta sqrt against tiny negative floats: sqrt(abs(...)) - Unify GMRES Hessenberg h_np assignment to avoid .real stripping producing wrong dtype for complex systems - Fix float() cast on dpnp scalar norm inside GMRES inner h_j1 line

…failures) The committed code used hypot(gbar, oldb) as delta_k which is the gamma (norm) from the PREVIOUS rotation step, not the correct diagonal entry from applying the previous Givens rotation to the current column. The correct Paige-Saunders (1975) two-rotation recurrence is: oldeps = epsln delta = cs * dbar + sn * alpha # apply previous rotation gbar_k = sn * dbar - cs * alpha # residual -> new rotation input epsln = sn * beta dbar = -cs * beta gamma = hypot(gbar_k, beta) # NEW rotation eliminates beta cs = gbar_k / gamma sn = beta / gamma w_new = (v - oldeps*w - delta*w2) / gamma # three-term update This matches scipy.sparse.linalg.minres and Choi (2006) eq. 6.11. The buggy recurrence produced solutions ~1.08x away from the true solution (rel_err ~1e0) instead of the expected ~1e-13. Co-authored-by: fix-minres-recurrence

…agusetty/dpnp into feature-sparse-linalg-solvers

Co-authored-by: Anton <100830759+antonwolfy@users.noreply.github.com>

…agusetty/dpnp into feature-sparse-linalg-solvers

This reverts commit 4074478.

antonwolfy · 2026-07-07T10:15:57Z

+and a small set of Krylov solvers (``cg``, ``gmres``, ``minres``).
+"""
+
+from __future__ import annotations


Do we need that? Is it used somewhere?

It is removed!

antonwolfy · 2026-07-07T10:31:41Z

+
+__all__ = [
+    "linalg",
+    "SparseABC",


Why do we need to expose SparseABC as part of scipy.sparse submodule?
SciPy does not did that.

antonwolfy · 2026-07-07T10:34:21Z

@@ -0,0 +1,27 @@
+"""Sparse base class and predicate, mirroring scipy/_lib/_sparse.py.


Should we mimic and move the file to dpnp/scipy/_lib/_sparse.py?

Moved it to mimic scipy style. Good idea indeed

antonwolfy · 2026-07-07T10:35:14Z

+
+from abc import ABC
+
+


Suggested change

__all__ = ["SparseABC", "issparse"]

antonwolfy · 2026-07-07T10:41:02Z

+    """Abstract base for all dpnp.scipy.sparse format classes."""
+
+
+def issparse(x):


We need to add rendering of documentation for the new sparse submodule.
Use doc/reference/scipy_linalg.rst as a reference.

antonwolfy · 2026-07-07T12:01:14Z

+        # set_csr_data + optimize_gemv must complete before any compute
+        # call can dispatch against the handle. This is the only blocking
+        # sync; subsequent matvecs return without waiting.
+        ev.wait()


Why do we need to have that call as a blocking? We can add the returned event under control of SequentialOrderManager.
Or alternatively we no need to return the event from _sparse_gemv_init and to implement the wait inside.

Added a aadd_event_pair to be consistent with the rest of codebase. Remove the blocking sync

antonwolfy · 2026-07-07T12:07:40Z

+                _manager.add_event_pair(ht_ev, comp_ev)
+                return y
+
+        # Dense fallback. Materialises ``self`` once -- this path is


I wonder if it is a fare option to implement the fallback silently. It does not seem present in CuPy. CuPy raises the explicit exception for unsupported dtypes. While SciPy support any dtypes, since implements own sparse kernels.

In fallback case we have to allocate the full dense array, which may cause protentional OOM issue.
And so leaving the fallback implementation under the user control might be more friendly.

Also, as an option, we can implement upcasting of input arrays to support more wide range of dtypes.

Dense fallback removed and dot raises TypeError for unsupported value/index dtype (CuPy-style, no O(M·N) densify/OOM). A few tests were added too

antonwolfy · 2026-07-07T12:11:24Z

+            return
+
+        release_fn = getattr(si, "_sparse_gemv_release", None)
+        if release_fn is None:


When there might be no _sparse_gemv_release available?

antonwolfy · 2026-07-07T12:12:56Z

+
+        try:
+            # pylint: disable-next=not-callable
+            release_fn(self._spmv_exec_q, handle, [])


Why no dependencies need to be passed?

Also release_fn returns the SYCL event but seems no handling

…agusetty/dpnp into feature-sparse-linalg-solvers

Co-authored-by: Anton <100830759+antonwolfy@users.noreply.github.com>

- _base.py: add copyright header, __all__, scipy-style issparse docstring - __init__.py: drop SparseABC from public __all__ (scipy does not expose it) - _csr.py: add copyright header - _csr.py: accept dpnp.ndarray and usm_ndarray via is_supported_array_type - _csr.py: use issparse for the csr-copy construction branch - _csr.py: validate components with check_supported_arrays_type - _csr.py: check queues via get_execution_queue + ExecutionPlacementError - _csr.py: support empty csr_matrix((M, N)) and CSR shape inference - _csr.py: allocate output with empty_like(shape=) - _csr.py: chain init event via SequentialOrderManager instead of blocking wait - _csr.py: remove silent dense fallback; raise for unsupported dtype (CuPy-style) - _csr.py: handle the returned release event in __del__ and clarify guards - _csr.py: raise NotImplementedError for unsupported scipy/cupy CSR ops - linalg/__init__.py: drop unused __future__ annotations import - doc: add scipy_sparse reference page and wire into toctree - tests: cover usm_ndarray construction/dot, empty/inferred shape, unsupported dtype

- _csr.py: pass submitted_events as depends to _sparse_gemv_init - _csr.py: pass depends + register returned event on release in __del__ - _iterative.py: allocate matvec output via empty_like(shape=) - _iterative.py: drop stale _make_fast_matvec fallback comments

Mirror SciPy's scipy/_lib/_sparse.py layout for the shared SparseABC / issparse helpers. - add dpnp/scipy/_lib package - update imports in scipy.sparse __init__ and _csr - register dpnp.scipy._lib in setup.py

- _csr.py: asarray-normalize dense/component inputs to dpnp.ndarray so internal ops (sort_indices, toarray, dot) work uniformly - _csr.py, _iterative.py: pass dpnp.get_usm_ndarray(...) to the sparse pybind bindings, accepting both dpnp.ndarray and usm_ndarray

…mpute_hu

- _interface.py: LinearOperator.dot and aslinearoperator accept usm_ndarray via is_supported_array_type (keep numpy-reject message) - _interface.py: numpydoc for LinearOperator, matvec, rmatvec, matmat, rmatmat, dot, aslinearoperator (Parameters/Returns, scipy refs) - _iterative.py: scipy-style numpydoc for cg, gmres, minres - tests: cover aslinearoperator/matvec with usm_ndarray dense input

…agusetty/dpnp into feature-sparse-linalg-solvers

abagusetty and others added 30 commits April 2, 2026 12:52

Add sparse.linalg iterative solvers

a7c8d31

Add sparse.linalg iterative solvers

e3e9052

Add sparse.linalg iterative solvers

c0e9fb5

Add sparse.linalg iterative solvers

943c52f

Add sparse.linalg iterative solvers

1191f7e

Fix deprecated tol kwarg in SciPy host fallback

2384185

Add tests for scipy.sparse.linalg: LinearOperator, cg, gmres, minres

58cc44b

Fix implicit numpy conversion; use .asnumpy() for dpnp arrays

3a50062

add onemkl sparse gemv pybind logic

8c68d98

sparse: add pybind11 module, CMakeLists, and hook _sparse_impl into _…

0c4a888

…iterative.py

Remove redundant sparse_gemv.hpp passthrough header

4993120

sparse: gemv.cpp includes gemv.hpp directly

a7ddc1c

sparse: capture exec_q from CSR data at closure construction

14cb5c4

minor cleanup for sparse extensions

ed58333

Fix SyntaxError: remove stray backslash in aslinearoperator hasattr s…

0a32b57

…tring Line 477: `hasattr(A, "rmatmat\")` had a Markdown-escaped backslash leaked into the Python source, causing an unterminated string literal. Fixed to `hasattr(A, "rmatmat")`.

sparse/linalg: pure-GPU CG/GMRES/MINRES, drop all CPU fallback paths,…

4292518

… port SciPy corner cases

Fix dtype.char AttributeError on dpnp dtype objects in CG/GMRES/MINRES

125dab5

Fix M guard, MINRES betacheck decay and gbar init in Paige-Saunders r…

2d753cf

…ecurrence

abagusetty and others added 4 commits July 4, 2026 11:08

Added a few more tests based on the reviews

ec1bb1f

Merge branch 'master' into feature-sparse-linalg-solvers

b7dc740

few more comments

eaf7c2b

Merge branch 'feature-sparse-linalg-solvers' of https://github.com/ab…

d8e2940

…agusetty/dpnp into feature-sparse-linalg-solvers

antonwolfy reviewed Jul 6, 2026

View reviewed changes

Comment thread dpnp/backend/extensions/sparse/CMakeLists.txt Outdated

Update dpnp/backend/extensions/blas/gemv.cpp

90f2228

Co-authored-by: Anton <100830759+antonwolfy@users.noreply.github.com>

antonwolfy reviewed Jul 6, 2026

View reviewed changes

Comment thread dpnp/backend/extensions/sparse/gemv.cpp Outdated

abagusetty and others added 8 commits July 6, 2026 15:45

address the comments

04cdba6

Merge branch 'master' into feature-sparse-linalg-solvers

e8b301b

fix build

a0ae7c2

Merge branch 'feature-sparse-linalg-solvers' of https://github.com/ab…

1ace5d2

…agusetty/dpnp into feature-sparse-linalg-solvers

add sum_duplicates for csr

a44c218

minor cleanup

e367f9e

format files

4074478

Revert "format files"

598e11b

This reverts commit 4074478.

antonwolfy reviewed Jul 7, 2026

View reviewed changes

abagusetty and others added 13 commits July 7, 2026 07:25

fix formatting issues

1c7d63d

Merge branch 'feature-sparse-linalg-solvers' of https://github.com/ab…

bcf05be

…agusetty/dpnp into feature-sparse-linalg-solvers

Merge branch 'master' into feature-sparse-linalg-solvers

b824866

Update dpnp/scipy/sparse/_csr.py

ed3ecf6

Co-authored-by: Anton <100830759+antonwolfy@users.noreply.github.com>

sparse: move _base.py to dpnp/scipy/_lib/_sparse.py

fd1c27e

Mirror SciPy's scipy/_lib/_sparse.py layout for the shared SparseABC / issparse helpers. - add dpnp/scipy/_lib package - update imports in scipy.sparse __init__ and _csr - register dpnp.scipy._lib in setup.py

sparse: trim _lib/_sparse.py docstring to match scipy style

f6c867c

sparse: use get_execution_queue + ExecutionPlacementError in _make_co…

a5a5094

…mpute_hu

sparse: import dpnp.tensor (not dpctl.tensor) to fix ModuleNotFoundError

e4647e3

Merge branch 'feature-sparse-linalg-solvers' of https://github.com/ab…

3762a9c

…agusetty/dpnp into feature-sparse-linalg-solvers

		@@ -0,0 +1,27 @@
		"""Sparse base class and predicate, mirroring scipy/_lib/_sparse.py.

		"""Abstract base for all dpnp.scipy.sparse format classes."""


		def issparse(x):

Uh oh!

Conversation

abagusetty commented Apr 9, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

antonwolfy Jul 7, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

5 participants

abagusetty commented Apr 9, 2026 •

edited

Loading

antonwolfy Jul 7, 2026 •

edited

Loading