Feature sparse linalg solvers

[abagusetty]

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?

[abagusetty]

@antonwolfy Apologies for the delay in getting back to your comments. All of them were addressed, whenever you get a chance. Thanks again!

[antonwolfy]

It should be a default floating point data type for the device on which the input data is allocated

[antonwolfy]

Missing description of attributes

[antonwolfy]

Suggested change

	functionality on top of DPNP arrays.
	functionality on top of DPNP arrays.

[antonwolfy]

It seems that has to be moved after the attributes block, otherwise rendered badly

[antonwolfy]

to ensure dtype is supported:

Suggested change

	dtype = dpnp.dtype(dtype)
	dtype = dpnp.empty(0, dtype=dtype).dtype

[antonwolfy]

dpnp might pass 7 to oneMKL -> OOB read of length-3 x for below example:

sp.csr_matrix((np.array([1.,2.]), np.array([0,7],np.int32), np.array([0,1,2],np.int32)), shape=(2,3))
# Out: <2x3 csr_matrix of dtype float64 with 2 stored elements>

while scipy:

a = sp.csr_matrix((np.array([1.,2.]), np.array([0,7],np.int32), np.array([0,1,2],np.int32)), shape=(2,3))
a.check_format()
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In[9], line 1
----> 1 a.check_format()

File /localdisk/work/antonvol/miniforge3/envs/dpnp_dev/lib/python3.13/site-packages/scipy/sparse/_compressed.py:212, in _cs_matrix.check_format(self, full_check)
    210 if self.nnz > 0:
    211     if self.indices.max() >= N:
--> 212         raise ValueError(f"indices must be < {N}")
    213     if self.indices.min() < 0:
    214         raise ValueError("indices must be >= 0")

ValueError: indices must be < 3

[antonwolfy]

• scipy: silent at construction, ValueError in check_format
• dpnp: row_lengths = diff = [2,-1] -> negative count to dpnp.repeat / bad offsets to oneMKL

sp.csr_matrix((np.array([1.,2.,3.]), np.array([0,1,2],np.int32), np.array([0,2,1],np.int32)), shape=(2,3))

[abagusetty]

Fixed it to be similar like scipy. It is silent at construction but check_format runs the checks lazily on the first use now

[antonwolfy]

alpha/beta passed by address of stack, because oneMath spmv takes const void* alpha and so might be destroyed when gemv_compute_impl returns

[abagusetty]

Fixed all the instance of this!

[antonwolfy]

self.A.dot(X[:, i]) on a C-contiguous 2-D X passes a non-unit-stride view.
sparse_gemv_compute only checks ndim/shape/dtype and reads x.get_data() contiguously → wrong results for L @ X/L.matmat(X) with ≥2 columns, no error.

[abagusetty]

Fixed. Thanks for catching this. sparse_gemv_compute now rejects non-contiguous x/y rather than silently consuming it. Added few unit-tests to cover this!

[antonwolfy]

The below example raises exception here:

import dpnp as np
import dpnp.scipy.sparse as sp

b0 = np.zeros(4, dtype=np.float64)
A_empty = sp.csr_matrix((4, 4), dtype=np.float64)

sp.linalg.cg(A_empty, b0, rtol=1e-8)
---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)
Cell In[12], line 1
----> 1 sp.linalg.cg(A_empty, b0, rtol=1e-8)

File /localdisk/work/antonvol/code/dpnp/dpnp/scipy/sparse/linalg/_iterative.py:373, in cg(A, b, x0, rtol, tol, maxiter, M, callback, atol)
    365     warnings.warn(
    366         "'tol' is deprecated in favor of 'rtol' and will be removed in "
    367         "a future release; use 'rtol' instead.",
    368         DeprecationWarning,
    369         stacklevel=2,
    370     )
    371     rtol = tol
--> 373 A_op, M_op, x, b, dtype = _make_system(A, M, x0, b)
    374 n = b.shape[0]
    376 bnrm = dpnp.linalg.norm(b)

File /localdisk/work/antonvol/code/dpnp/dpnp/scipy/sparse/linalg/_iterative.py:271, in _make_system(A, M, x0, b)
    268         M_op = _FastMOp()
    270 # Inject fast CSR SpMV for A if available.
--> 271 fast_mv = _make_fast_matvec(A)
    272 if fast_mv is not None:
    273     _orig = A_op

File /localdisk/work/antonvol/code/dpnp/dpnp/scipy/sparse/linalg/_iterative.py:178, in _make_fast_matvec(A)
    176     return None
    177 # pylint: disable-next=protected-access
--> 178 handle_info = A._ensure_spmv_handle()
    179 if handle_info is None:
    180     return None

File /localdisk/work/antonvol/code/dpnp/dpnp/scipy/sparse/_csr.py:445, in csr_matrix._ensure_spmv_handle(self)
    443 _manager = _dpu.SequentialOrderManager[exec_q]
    444 # pylint: disable-next=protected-access
--> 445 handle, val_type_id, ev = _si._sparse_gemv_init(
    446     exec_q,
    447     0,  # trans=N (forward)
    448     _dpnp.get_usm_ndarray(self.indptr),
    449     _dpnp.get_usm_ndarray(self.indices),
    450     _dpnp.get_usm_ndarray(self.data),
    451     int(self._shape[0]),
    452     int(self._shape[1]),
    453     int(self.data.shape[0]),
    454     _manager.submitted_events,
    455 )
    457 # set_csr_data + optimize_gemv must complete before the first
    458 # compute; chain the init event through the queue's order manager
    459 # so the first _sparse_gemv_compute depends on it (non-blocking).
    460 _manager.add_event_pair(ev, ev)

RuntimeError: sparse_gemv_init: MKL exception in set_csr_data: oneapi::mkl::sparse::set_csr_data: invalid argument: nnz is = 0 provided with at least one of `ja`/`a` != nullptr

but passes with scipy

[abagusetty]

Fixed!

[abagusetty]

@antonwolfy Is there a timeline for the 0.21.0 release/code freeze. Just wanted to prioritize this PR accordingly on my end

[antonwolfy]

@antonwolfy Is there a timeline for the 0.21.0 release/code freeze. Just wanted to prioritize this PR accordingly on my end

@abagusetty, as for now, we are planning the freeze by mid of September.

[antonwolfy]

That causes tests failure on device without fp64:

FAILED tests/test_scipy_sparse_linalg.py::TestLinearOperatorAlgebra::test_scaled_linear_operator[float32] - ValueError: Device Intel(R) Iris(R) Xe Graphics does not provide native support for double-precision floating point type.
--
tests/test_scipy_sparse_linalg.py::TestLinearOperatorAlgebra::test_scaled_linear_operator[complex64] - ValueError: Device Intel(R) Iris(R) Xe Graphics does not provide native support for double-precision floating point type.

[vlad-perevezentsev]

It looks like non-contiguous inputs are not supported here.
sparse_gemv_compute reads x and y linearly but does not check that they are contiguous.
The same applies to row_ptr, col_ind and values in sparse_gemv_init

For strided inputs this can silently produce incorrect results instead of raising an error

In [1]: import dpnp

In [2]: A = dpnp.scipy.sparse.csr_matrix(dpnp.array([[1., 2., 0.], [0., 3., 4.], [5., 0., 6.]]))

In [3]: x = dpnp.arange(6, dtype='f8')[::2]

In [4]: A.dot(x)
Out[4]: array([ 2., 11., 12.])

# scipy 
In [1]: import scipy

In [2]: import numpy

In [3]: A = scipy.sparse.csc_matrix(numpy.array([[1., 2., 0.], [0., 3., 4.], [5., 0., 6.]]))

In [4]: x = numpy.arange(6, dtype=numpy.float64)[::2]

In [5]: A.dot(x)
Out[5]: array([ 4., 22., 24.])

Also the current tests do not cover non-contiguous inputs

[vlad-perevezentsev]

Would it make sense to make them read-only properties?
The oneMKL handle is cached using pointers to these buffers but reassigning any of them does not invalidate the handle and may leave it pointing to freed USM memory

In [1]: import numpy

In [2]: import dpnp

In [3]: A = dpnp.scipy.sparse.csr_matrix(dpnp.array(numpy.diag([1., 2., 3.])))

In [4]: A.dot(dpnp.ones(3))
Out[4]: array([1., 2., 3.])

In [5]: A.data = dpnp.array([10.,20.,30.])

In [6]: A.dot(dpnp.ones(3))
Out[6]: array([1., 1., 1.])

#scipy

In [1]: import scipy

In [2]: import numpy

In [3]: A = scipy.sparse.csr_matrix(numpy.diag([1., 2., 3.]))

In [4]: A.dot(numpy.ones(3))
Out[4]: array([1., 2., 3.])

In [5]: A.data = numpy.array([10.,20.,30.])

In [6]: A.dot(numpy.ones(3))
Out[6]: array([10., 20., 30.])

[vlad-perevezentsev]

float64 and complex128 may not be supported on GPUs without fp64

[vlad-perevezentsev]

oneMKL set_csr_data rejects nnz == 0 line 506 in _csr.py
Should we add the same guard here?

[vlad-perevezentsev]

This check should be moved above shape = tuple(int(s) for s in shape)

Right now (3.7, 3.2) is converted to (3, 3) before _isshape() runs,so non-integer dimensions are silently accepted

In [5]: from dpnp.scipy.sparse.linalg import LinearOperator

In [6]: LinearOperator((3.7, 3.2), matvec=lambda v: v).shape
Out[6]: (3, 3)

# scipy
In [10]: from scipy.sparse.linalg import LinearOperator

In [11]: LinearOperator((3.7, 3.2), matvec=lambda v: v).shape

ValueError: invalid shape (3.7, 3.2)

[vlad-perevezentsev]

Suggested change

	return b, 0
	return dpnp.zeros_like(b), 0

Just like in def cg

[vlad-perevezentsev]

Would not it be simpler to just use self.A.dot(x) here for both branches?
dpnp.ndarray.dot with 1-D argument already dispatches to gemv