Build with sycl_thrust for device examples if building oneMKL.

CMakeLists.txt

@@ -20,9 +20,16 @@ endif()
 include(FetchContent)
 
 if (ENABLE_ONEMKL_SYCL)
+  set(SPBLAS_GPU_BACKEND ON)
   find_package(MKL REQUIRED)
   target_link_libraries(spblas INTERFACE MKL::MKL_SYCL) # SYCL APIs
   set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -DSPBLAS_ENABLE_ONEMKL_SYCL")
+
+  FetchContent_Declare(
+  sycl_thrust
+  GIT_REPOSITORY https://github.com/SparseBLAS/sycl-thrust.git
+  GIT_TAG main)
+  FetchContent_MakeAvailable(sycl_thrust)
 endif()
 
 if (ENABLE_ARMPL)

examples/CMakeLists.txt

@@ -11,7 +11,6 @@ if (NOT SPBLAS_GPU_BACKEND)
   add_example(matrix_opt_example)
   add_example(spmm_csc)
 else()
-  find_package(rocthrust REQUIRED)
   add_subdirectory(device)
 endif()
 

examples/device/CMakeLists.txt

@@ -1,12 +1,18 @@
 function(add_device_example example_name)
+  add_executable(${example_name} ${example_name}.cpp)
+
   if (ENABLE_ROCSPARSE)
+    find_package(rocthrust REQUIRED)
     set_source_files_properties(${example_name}.cpp PROPERTIES LANGUAGE HIP)
+    target_link_libraries(${example_name} rocthrust)
   # elseif (ENABLE_CUSPARSE)
     # cuSPARSE linking details will go here.
+  elseif (ENABLE_ONEMKL_SYCL)
+    target_link_libraries(${example_name} sycl_thrust)
   else()
     message(FATAL_ERROR "Device backend not found.")
   endif()
-  add_executable(${example_name} ${example_name}.cpp)
+
   target_link_libraries(${example_name} spblas fmt)
 endfunction()
 

examples/device/simple_spmv.cpp

@@ -56,8 +56,7 @@ int main(int argc, char** argv) {
   std::span<value_t> y_span(d_y.data().get(), m);
 
   // y = A * x
-  spblas::spmv_state_t state;
-  spblas::multiply(state, a, x_span, y_span);
+  spblas::multiply(a, x_span, y_span);
 
   thrust::copy(d_y.begin(), d_y.end(), y.begin());
 

include/spblas/vendor/onemkl_sycl/detail/detail.hpp

@@ -1,4 +1,6 @@
 #pragma once
 
 #include "create_matrix_handle.hpp"
+#include "execution_policy.hpp"
 #include "get_matrix_handle.hpp"
+#include "get_queue.hpp"

include/spblas/vendor/onemkl_sycl/detail/execution_policy.hpp

@@ -0,0 +1,52 @@
+#pragma once
+
+#include <spblas/vendor/onemkl_sycl/detail/get_pointer_device.hpp>
+#include <sycl/sycl.hpp>
+
+namespace spblas {
+
+namespace mkl {
+
+class parallel_policy {
+public:
+  parallel_policy() {}
+
+  template <typename T>
+  sycl::queue get_queue(T* ptr) const {
+    return spblas::__mkl::get_pointer_queue(ptr);
+  }
+
+  sycl::queue get_queue() const {
+    return sycl::queue(sycl::default_selector_v);
+  }
+};
+
+class device_policy {
+public:
+  device_policy(const sycl::queue& queue) : queue_(queue) {}
+
+  sycl::queue& get_queue() {
+    return queue_;
+  }
+
+  const sycl::queue& get_queue() const {
+    return queue_;
+  }
+
+  sycl::device get_device() const {
+    return queue_.get_device();
+  }
+
+  sycl::context get_context() const {
+    return queue_.get_context();
+  }
+
+private:
+  sycl::queue queue_;
+};
+
+inline parallel_policy par;
+
+} // namespace mkl
+
+} // namespace spblas

include/spblas/vendor/onemkl_sycl/detail/get_pointer_device.hpp

@@ -0,0 +1,46 @@
+#pragma once
+
+#include <sycl/sycl.hpp>
+#include <vector>
+
+namespace spblas {
+
+namespace __mkl {
+
+inline std::vector<sycl::context> global_contexts_;
+
+template <typename T>
+std::pair<sycl::device, sycl::context> get_pointer_device(T* ptr) {
+  if (global_contexts_.empty()) {
+    for (auto&& platform : sycl::platform::get_platforms()) {
+      sycl::context context(platform.get_devices());
+
+      global_contexts_.push_back(context);
+    }
+  }
+
+  for (auto&& context : global_contexts_) {
+    try {
+      sycl::device device = sycl::get_pointer_device(ptr, context);
+      return {device, context};
+    } catch (...) {
+    }
+  }
+
+  throw std::runtime_error(
+      "get_pointer_device: could not locate device corresponding to pointer");
+}
+
+template <typename T>
+sycl::queue get_pointer_queue(T* ptr) {
+  try {
+    auto&& [device, context] = get_pointer_device(ptr);
+    return sycl::queue(context, device);
+  } catch (...) {
+    return sycl::queue(sycl::cpu_selector_v);
+  }
+}
+
+} // namespace __mkl
+
+} // namespace spblas

include/spblas/vendor/onemkl_sycl/detail/get_queue.hpp

@@ -0,0 +1,40 @@
+#pragma once
+
+#include <spblas/vendor/onemkl_sycl/detail/execution_policy.hpp>
+
+namespace spblas {
+
+namespace __mkl {
+
+template <typename T>
+sycl::queue get_queue(const spblas::mkl::parallel_policy& policy, T* ptr) {
+  return policy.get_queue(ptr);
+}
+
+template <typename T>
+sycl::queue& get_queue(spblas::mkl::device_policy& policy, T* ptr) {
+  return policy.get_queue();
+}
+
+} // namespace __mkl
+
+} // namespace spblas
+
+#if __has_include(<thrust/execution_policy.h>)
+
+#include <thrust/execution_policy.h>
+
+namespace spblas {
+
+namespace __mkl {
+
+template <typename T>
+sycl::queue& get_queue(thrust::execution_policy& policy, T* ptr) {
+  return policy.get_queue();
+}
+
+} // namespace __mkl
+
+} // namespace spblas
+
+#endif

include/spblas/vendor/onemkl_sycl/spmv_impl.hpp

@@ -26,18 +26,20 @@
 
 namespace spblas {
 
-template <matrix A, vector X, vector Y>
+template <typename ExecutionPolicy, matrix A, vector X, vector Y>
   requires((__detail::has_csr_base<A> || __detail::has_csc_base<A>) &&
            __detail::has_contiguous_range_base<X> &&
            __ranges::contiguous_range<Y>)
-void multiply(A&& a, X&& x, Y&& y) {
+void multiply(ExecutionPolicy&& policy, A&& a, X&& x, Y&& y) {
   log_trace("");
   auto x_base = __detail::get_ultimate_base(x);
 
   auto alpha_optional = __detail::get_scaling_factor(a, x);
   tensor_scalar_t<A> alpha = alpha_optional.value_or(1);
 
-  sycl::queue q(sycl::cpu_selector_v);
+  auto a_data = __detail::get_ultimate_base(a).values().data();
+
+  auto&& q = __mkl::get_queue(policy, a_data);
 
   auto a_handle = __mkl::get_matrix_handle(q, a);
   auto a_transpose = __mkl::get_transpose(a);
@@ -51,4 +53,13 @@ void multiply(A&& a, X&& x, Y&& y) {
   }
 }
 
+template <matrix A, vector X, vector Y>
+  requires((__detail::has_csr_base<A> || __detail::has_csc_base<A>) &&
+           __detail::has_contiguous_range_base<X> &&
+           __ranges::contiguous_range<Y>)
+void multiply(A&& a, X&& x, Y&& y) {
+  multiply(mkl::par, std::forward<A>(a), std::forward<X>(x),
+           std::forward<Y>(y));
+}
+
 } // namespace spblas

test/gtest/CMakeLists.txt

@@ -1,23 +1,39 @@
 enable_testing()
 
+set(TEST_SOURCES)
+
 if (NOT ENABLE_ROCSPARSE)
-  add_executable(
-    spblas-tests
-    spmv_test.cpp
-    spmm_test.cpp
-    spgemm_test.cpp
-    spgemm_csr_csc.cpp
-    add_test.cpp
-    transpose_test.cpp
-    triangular_solve_test.cpp
-  )
-elseif(ENABLE_ROCSPARSE)
-  set_source_files_properties(rocsparse/spmv_test.cpp PROPERTIES LANGUAGE HIP)
-  add_executable(spblas-tests
-                 rocsparse/spmv_test.cpp)
+  list(APPEND TEST_SOURCES
+       spmv_test.cpp
+       spmm_test.cpp
+       spgemm_test.cpp
+       spgemm_csr_csc.cpp
+       add_test.cpp
+       transpose_test.cpp
+       triangular_solve_test.cpp)
 endif()
 
+if (SPBLAS_GPU_BACKEND)
+  list(APPEND TEST_SOURCES
+       thrust/spmv_test.cpp)
+
+  if (ENABLE_ROCSPARSE)
+    set_source_files_properties(${TEST_SOURCES} PROPERTIES LANGUAGE HIP)
+  endif()
+endif()
+
+add_executable(spblas-tests ${TEST_SOURCES})
+
 target_link_libraries(spblas-tests spblas fmt GTest::gtest_main)
 
+if (ENABLE_ROCSPARSE)
+  find_package(rocthrust REQUIRED)
+  target_link_libraries(spblas-tests rocthrust)
+endif()
+
+if (ENABLE_ONEMKL_SYCL)
+  target_link_libraries(spblas-tests sycl_thrust)
+endif()
+
 include(GoogleTest)
 gtest_discover_tests(spblas-tests)

test/gtest/rocsparse/spmv_test.cpp → test/gtest/thrust/spmv_test.cpp

@@ -1,4 +1,3 @@
-
 #include "../util.hpp"
 #include <spblas/spblas.hpp>
 
@@ -9,7 +8,7 @@ using value_t = float;
 using index_t = spblas::index_t;
 using offset_t = spblas::offset_t;
 
-TEST(CsrView, SpMV) {
+TEST(thrust_CsrView, SpMV) {
   for (auto&& [num_rows, num_cols, nnz] : util::dims) {
     auto [values, rowptr, colind, shape, _] =
         spblas::generate_csr<value_t, index_t, offset_t>(num_rows, num_cols,
@@ -32,8 +31,7 @@ TEST(CsrView, SpMV) {
     std::span<value_t> b_span(d_b.data().get(), num_cols);
     std::span<value_t> c_span(d_c.data().get(), num_rows);
 
-    spblas::spmv_state_t state;
-    spblas::multiply(state, a, b_span, c_span);
+    spblas::multiply(a, b_span, c_span);
 
     thrust::copy(d_c.begin(), d_c.end(), c.begin());
 
@@ -53,7 +51,7 @@ TEST(CsrView, SpMV) {
   }
 }
 
-TEST(CsrView, SpMV_Ascaled) {
+TEST(thrust_CsrView, SpMV_Ascaled) {
   for (auto&& [num_rows, num_cols, nnz] :
        {std::tuple(1000, 100, 100), std::tuple(100, 1000, 10000),
         std::tuple(40, 40, 1000)}) {
@@ -79,8 +77,7 @@ TEST(CsrView, SpMV_Ascaled) {
       std::span<value_t> b_span(d_b.data().get(), num_cols);
       std::span<value_t> c_span(d_c.data().get(), num_rows);
 
-      spblas::spmv_state_t state;
-      spblas::multiply(state, spblas::scaled(alpha, a), b_span, c_span);
+      spblas::multiply(spblas::scaled(alpha, a), b_span, c_span);
 
       thrust::copy(d_c.begin(), d_c.end(), c.begin());
 
@@ -101,7 +98,7 @@ TEST(CsrView, SpMV_Ascaled) {
   }
 }
 
-TEST(CsrView, SpMV_BScaled) {
+TEST(thrust_CsrView, SpMV_BScaled) {
   for (auto&& [num_rows, num_cols, nnz] :
        {std::tuple(1000, 100, 100), std::tuple(100, 1000, 10000),
         std::tuple(40, 40, 1000)}) {
@@ -127,8 +124,7 @@ TEST(CsrView, SpMV_BScaled) {
       std::span<value_t> b_span(d_b.data().get(), num_cols);
       std::span<value_t> c_span(d_c.data().get(), num_rows);
 
-      spblas::spmv_state_t state;
-      spblas::multiply(state, a, spblas::scaled(alpha, b_span), c_span);
+      spblas::multiply(a, spblas::scaled(alpha, b_span), c_span);
 
       thrust::copy(d_c.begin(), d_c.end(), c.begin());