Feature sparse linalg solvers

CHANGELOG.md

@@ -13,6 +13,7 @@ This release is compatible with NumPy 2.4.5.
 * Added C API functions for `dpnp.tensor.usm_ndarray` setters and getters to avoid ABI breakage if `dpnp.tensor.usm_ndarray` is modified [gh-2866](https://github.com/IntelPython/dpnp/pull/2866)
 * Added support for buffer protocol objects as advanced index keys in `dpnp.ndarray` [#2889](https://github.com/IntelPython/dpnp/pull/2889)
 * Added `--includes` and `--include-dir` options to the `dpnp` CLI [#2916](https://github.com/IntelPython/dpnp/pull/2916)
+* Added implementation of `dpnp.scipy.sparse.linalg import LinearOperator, cg, gmres, minres` [#2841](https://github.com/IntelPython/dpnp/pull/2841)
 
 ### Changed
 

conda-recipe/meta.yaml

@@ -50,6 +50,7 @@ requirements:
       - {{ pin_compatible('onemkl-sycl-lapack', min_pin='x.x', max_pin='x') }}
       - {{ pin_compatible('onemkl-sycl-rng', min_pin='x.x', max_pin='x') }}
       - {{ pin_compatible('onemkl-sycl-vm', min_pin='x.x', max_pin='x') }}
+      - {{ pin_compatible('onemkl-sycl-sparse', min_pin='x.x', max_pin='x') }}
       - numpy
       - intel-gpu-ocl-icd-system
 

dpnp/CMakeLists.txt

@@ -185,6 +185,7 @@ add_subdirectory(backend/extensions/statistics)
 add_subdirectory(backend/extensions/ufunc)
 add_subdirectory(backend/extensions/vm)
 add_subdirectory(backend/extensions/window)
+add_subdirectory(backend/extensions/sparse)
 
 add_subdirectory(dpnp_algo)
 add_subdirectory(dpnp_utils)

dpnp/backend/extensions/blas/blas_py.cpp

@@ -151,6 +151,32 @@ PYBIND11_MODULE(_blas_impl, m)
               py::arg("depends") = py::list());
     }
 
+    {
+        // y = alpha * op(A) * x + beta * y, the full BLAS gemv form.
+        // Used by dpnp.scipy.sparse.linalg.gmres to fuse the Arnoldi
+        // orthogonalisation (u -= V @ h) into a single oneMKL call
+        // (trans_op=0, alpha=-1, beta=1) and to write the Hessenberg
+        // column h = V^H @ u directly into the Hessenberg matrix slice
+        // (trans_op=2, alpha=1, beta=0).
+        //
+        // trans_op selects the op applied to A:
+        //   0 = N   y = alpha * A   @ x + beta * y
+        //   1 = T   y = alpha * A^T @ x + beta * y
+        //   2 = C   y = alpha * A^H @ x + beta * y   (F-contig only)
+        //
+        // For complex matrices the scalars must be exactly
+        // representable as their real form: callers pass {-1, 0, 1};
+        // fractional or complex scalars would lose the imaginary
+        // component on the C++ cast.
+        m.def("_gemv_alpha_beta", &blas_ns::gemv_alpha_beta,
+              "Call `gemv` from oneMKL BLAS library with explicit "
+              "alpha, beta, and tri-state trans_op (0=N, 1=T, 2=C). "
+              "Computes y = alpha * op(A) * x + beta * y.",
+              py::arg("sycl_queue"), py::arg("matrixA"), py::arg("vectorX"),
+              py::arg("vectorY"), py::arg("trans_op"), py::arg("alpha"),
+              py::arg("beta"), py::arg("depends") = py::list());
+    }
+
     {
         m.def("_syrk", &blas_ns::syrk,
               "Call `syrk` from oneMKL BLAS library to compute "

dpnp/backend/extensions/blas/gemv.cpp

@@ -49,14 +49,20 @@ namespace type_utils = dpnp::tensor::type_utils;
 
 using ext::common::init_dispatch_vector;
 
+// Impl signature now carries alpha and beta as double. Each per-T impl
+// casts them to T at the very end (T(alpha_d), T(beta_d)) so the same
+// dispatch vector serves both the legacy alpha=1/beta=0 wrapper and
+// the new alpha_beta entry point used by the GMRES Arnoldi fast path.
 typedef sycl::event (*gemv_impl_fn_ptr_t)(sycl::queue &,
                                           oneapi::mkl::transpose,
                                           const std::int64_t,
                                           const std::int64_t,
+                                          const double, // alpha
                                           const char *,
                                           const std::int64_t,
                                           const char *,
                                           const std::int64_t,
+                                          const double, // beta
                                           char *,
                                           const std::int64_t,
                                           const bool,
@@ -69,10 +75,12 @@ static sycl::event gemv_impl(sycl::queue &exec_q,
                              oneapi::mkl::transpose transA,
                              const std::int64_t m,
                              const std::int64_t n,
+                             const double alpha_d,
                              const char *matrixA,
                              const std::int64_t lda,
                              const char *vectorX,
                              const std::int64_t incx,
+                             const double beta_d,
                              char *vectorY,
                              const std::int64_t incy,
                              const bool is_row_major,
@@ -84,6 +92,15 @@ static sycl::event gemv_impl(sycl::queue &exec_q,
     const T *x = reinterpret_cast<const T *>(vectorX);
     T *y = reinterpret_cast<T *>(vectorY);
 
+    // Cast alpha/beta into the matrix value type. For complex T the
+    // single-argument constructor sets the imaginary component to
+    // zero, which is exact for the GMRES use case (alpha and beta are
+    // always one of {-1, 0, 1}) and for the dpnp.dot wrapper (alpha=1,
+    // beta=0). Callers passing fractional or complex scalars through
+    // this path would lose the imaginary component silently.
+    const T alpha = static_cast<T>(alpha_d);
+    const T beta = static_cast<T>(beta_d);
+
     std::stringstream error_msg;
     bool is_exception_caught = false;
 
@@ -112,13 +129,13 @@ static sycl::event gemv_impl(sycl::queue &exec_q,
                     // or 'C' for a conjugate transpose.
             m,      // Number of rows in matrix A.
             n,      // Number of columns in matrix A.
-            T(1),   // Scaling factor for the matrix-vector product.
+            alpha,  // Scaling factor for the matrix-vector product.
             a,      // Pointer to the input matrix A.
             lda,    // Leading dimension of matrix A, which is the
                     // stride between successive rows (for row major layout).
             x,      // Pointer to the input vector x.
             incx,   // The stride of vector x.
-            T(0),   // Scaling factor for vector y.
+            beta,   // Scaling factor for vector y.
             y,      // Pointer to output vector y, where the result is stored.
             incy,   // The stride of vector y.
             depends);
@@ -141,14 +158,38 @@ static sycl::event gemv_impl(sycl::queue &exec_q,
     return gemv_event;
 }
 
-std::pair<sycl::event, sycl::event>
-    gemv(sycl::queue &exec_q,
-         const dpnp::tensor::usm_ndarray &matrixA,
-         const dpnp::tensor::usm_ndarray &vectorX,
-         const dpnp::tensor::usm_ndarray &vectorY,
-         const bool transpose,
-         const std::vector<sycl::event> &depends)
+// Shared validation + dispatch. Both gemv() (alpha=1, beta=0) and
+// gemv_alpha_beta() funnel through here so behaviour stays identical
+// across the two entry points.
+//
+// ``trans_op`` is a tri-state matching oneapi::mkl::transpose:
+//      0 = N (no transpose),
+//      1 = T (plain transpose),
+//      2 = C (conjugate-transpose, complex only).
+//
+// The legacy gemv() entry-point only ever needs N/T (real or complex,
+// no conjugate semantics) and forwards trans_op = 0 or 1. The
+// gemv_alpha_beta() entry-point exposes the full tri-state so the
+// GMRES Arnoldi inner step can request V^H directly instead of
+// post-conjugating the result of a T-mode gemv -- which is
+// mathematically wrong for a complex right-hand vector (the identity
+// conj(V^T @ u) == V^H @ u holds only when u is real-valued).
+static std::pair<sycl::event, sycl::event>
+    gemv_dispatch(sycl::queue &exec_q,
+                  const dpnp::tensor::usm_ndarray &matrixA,
+                  const dpnp::tensor::usm_ndarray &vectorX,
+                  const dpnp::tensor::usm_ndarray &vectorY,
+                  const int trans_op,
+                  const double alpha,
+                  const double beta,
+                  const std::vector<sycl::event> &depends)
 {
+    if (trans_op < 0 || trans_op > 2) {
+        throw py::value_error("gemv: trans_op must be 0 (N), 1 (T), or 2 (C).");
+    }
+    const bool is_transposed = (trans_op != 0);
+    const bool is_conj_trans = (trans_op == 2);
+
     const int matrixA_nd = matrixA.get_ndim();
     const int vectorX_nd = vectorX.get_ndim();
     const int vectorY_nd = vectorY.get_ndim();
@@ -182,10 +223,22 @@ std::pair<sycl::event, sycl::event>
             "Input matrix is not c-contiguous nor f-contiguous.");
     }
 
+    // Conjugate-transpose is only wired up for column-major (F-contig)
+    // matrices. The row-major remap (treating a C-contig matrix as its
+    // column-major transpose) does not extend cleanly to the C op
+    // because (A^T)^H == conj(A), which oneMKL does not expose as a
+    // gemv mode. Callers needing C-mode on row-major input must
+    // F-contigify first (e.g. via dpnp.asarray(A, order="F")).
+    if (is_conj_trans && !is_matrixA_f_contig) {
+        throw py::value_error(
+            "gemv: trans_op = 2 (conjugate-transpose) requires an "
+            "F-contiguous matrix; pass dpnp.asarray(A, order='F') first.");
+    }
+
     const py::ssize_t *a_shape = matrixA.get_shape_raw();
     const py::ssize_t *x_shape = vectorX.get_shape_raw();
     const py::ssize_t *y_shape = vectorY.get_shape_raw();
-    if (transpose) {
+    if (is_transposed) {
         if (a_shape[0] != x_shape[0]) {
             throw py::value_error("The number of rows in A must be equal to "
                                   "the number of elements in X.");
@@ -209,6 +262,9 @@ std::pair<sycl::event, sycl::event>
     oneapi::mkl::transpose transA;
     std::size_t src_nelems;
 
+    // Resolve the storage layout into the oneMKL transpose mode.
+    // Conjugate-transpose is constrained to F-contig above; the
+    // row-major branch therefore only sees N/T here.
 // cuBLAS supports only column-major storage
 #if defined(USE_ONEMATH_CUBLAS)
     constexpr bool is_row_major = false;
@@ -218,7 +274,11 @@ std::pair<sycl::event, sycl::event>
     if (is_matrixA_f_contig) {
         m = a_shape[0];
         n = a_shape[1];
-        if (transpose) {
+        if (is_conj_trans) {
+            transA = oneapi::mkl::transpose::C;
+            src_nelems = n;
+        }
+        else if (is_transposed) {
             transA = oneapi::mkl::transpose::T;
             src_nelems = n;
         }
@@ -228,9 +288,11 @@ std::pair<sycl::event, sycl::event>
         }
     }
     else {
+        // Row-major-as-column-major swap. is_conj_trans is rejected
+        // above, so only N/T need handling.
         m = a_shape[1];
         n = a_shape[0];
-        if (transpose) {
+        if (is_transposed) {
             transA = oneapi::mkl::transpose::N;
             src_nelems = m;
         }
@@ -248,7 +310,11 @@ std::pair<sycl::event, sycl::event>
     const std::int64_t m = a_shape[0];
     const std::int64_t n = a_shape[1];
 
-    if (transpose) {
+    if (is_conj_trans) {
+        transA = oneapi::mkl::transpose::C;
+        src_nelems = n;
+    }
+    else if (is_transposed) {
         transA = oneapi::mkl::transpose::T;
         src_nelems = n;
     }
@@ -299,16 +365,55 @@ std::pair<sycl::event, sycl::event>
         y_typeless_ptr -= (y_shape[0] - 1) * std::abs(incy) * y_elemsize;
     }
 
-    sycl::event gemv_ev =
-        gemv_fn(exec_q, transA, m, n, a_typeless_ptr, lda, x_typeless_ptr, incx,
-                y_typeless_ptr, incy, is_row_major, depends);
+    sycl::event gemv_ev = gemv_fn(exec_q, transA, m, n, alpha, a_typeless_ptr,
+                                  lda, x_typeless_ptr, incx, beta,
+                                  y_typeless_ptr, incy, is_row_major, depends);
 
     sycl::event args_ev = dpnp::utils::keep_args_alive(
         exec_q, {matrixA, vectorX, vectorY}, {gemv_ev});
 
     return std::make_pair(args_ev, gemv_ev);
 }
 
+std::pair<sycl::event, sycl::event>
+    gemv(sycl::queue &exec_q,
+         const dpnp::tensor::usm_ndarray &matrixA,
+         const dpnp::tensor::usm_ndarray &vectorX,
+         const dpnp::tensor::usm_ndarray &vectorY,
+         const bool transpose,
+         const std::vector<sycl::event> &depends)
+{
+    // Legacy alpha=1, beta=0 wrapper. Existing dpnp.dot callers expect
+    // this exact behaviour (N or plain T only, never conjugate-
+    // transpose), so we forward through the shared dispatch mapping
+    // the bool argument to the {0=N, 1=T} subset of the tri-state.
+    const int trans_op = transpose ? 1 : 0;
+    return gemv_dispatch(exec_q, matrixA, vectorX, vectorY, trans_op,
+                         /*alpha=*/1.0, /*beta=*/0.0, depends);
+}
+
+std::pair<sycl::event, sycl::event>
+    gemv_alpha_beta(sycl::queue &exec_q,
+                    const dpnp::tensor::usm_ndarray &matrixA,
+                    const dpnp::tensor::usm_ndarray &vectorX,
+                    const dpnp::tensor::usm_ndarray &vectorY,
+                    const int trans_op,
+                    const double alpha,
+                    const double beta,
+                    const std::vector<sycl::event> &depends)
+{
+    // Caller-supplied alpha / beta and full tri-state transpose.
+    // Used by the GMRES Arnoldi step to fuse  u -= V @ h
+    // (trans_op=0, alpha=-1, beta=1) into a single gemv kernel, and
+    // to write h = V^H @ u directly into a Hessenberg column slice
+    // (trans_op=2, alpha=1, beta=0). For complex matrices the scalars
+    // must be exactly representable as their real form -- callers
+    // pass 1/0/-1 only, see the impl-level comment for the imag-loss
+    // caveat.
+    return gemv_dispatch(exec_q, matrixA, vectorX, vectorY, trans_op, alpha,
+                         beta, depends);
+}
+
 template <typename fnT, typename varT>
 struct GemvContigFactory
 {

dpnp/backend/extensions/blas/gemv.hpp

@@ -35,6 +35,8 @@
 
 namespace dpnp::extensions::blas
 {
+// Convenience wrapper: alpha = 1, beta = 0. Preserved for the existing
+// dpnp call sites (dpnp.dot etc.) that do not care about scaling.
 extern std::pair<sycl::event, sycl::event>
     gemv(sycl::queue &exec_q,
          const dpnp::tensor::usm_ndarray &matrixA,
@@ -43,5 +45,36 @@ extern std::pair<sycl::event, sycl::event>
          const bool transpose,
          const std::vector<sycl::event> &depends);
 
+// Full y = alpha * op(A) * x + beta * y form. Required by the GMRES
+// Arnoldi step where we fuse u -= V @ h into a single gemv call
+// (alpha = -1, beta = 1), and where we write h = V^H @ u directly into
+// a Hessenberg column slice (alpha = 1, beta = 0). Both alpha and
+// beta arrive as double on the Python side and are cast to the matrix
+// value type inside the impl -- complex callers should use 1 / 0 / -1
+// (representable exactly) to avoid silent imaginary loss.
+//
+// ``trans_op`` selects the operation applied to A:
+//      0 = N  (no transpose)         y = alpha * A   @ x + beta * y
+//      1 = T  (transpose)            y = alpha * A^T @ x + beta * y
+//      2 = C  (conjugate-transpose)  y = alpha * A^H @ x + beta * y
+//
+// For real-valued A, T and C are equivalent. For complex A they
+// differ, and C is required for any algorithm that performs a
+// Hermitian inner product through gemv -- the GMRES Arnoldi step
+// (Gram-Schmidt over a complex Krylov basis) being the canonical
+// example. ``trans_op = 2`` is currently only supported for
+// F-contiguous (column-major) matrices; the row-major code path
+// for conjugate-transpose would require an explicit element-wise
+// conjugate pass and is not wired up here.
+extern std::pair<sycl::event, sycl::event>
+    gemv_alpha_beta(sycl::queue &exec_q,
+                    const dpnp::tensor::usm_ndarray &matrixA,
+                    const dpnp::tensor::usm_ndarray &vectorX,
+                    const dpnp::tensor::usm_ndarray &vectorY,
+                    const int trans_op,
+                    const double alpha,
+                    const double beta,
+                    const std::vector<sycl::event> &depends);
+
 extern void init_gemv_dispatch_vector(void);
 } // namespace dpnp::extensions::blas