Add sparse SYCL matmatd primitive

cpp/oneapi/dal/backend/primitives/sparse_blas.hpp

@@ -19,3 +19,4 @@
 #include "oneapi/dal/backend/primitives/sparse_blas/set_csr_data.hpp"
 #include "oneapi/dal/backend/primitives/sparse_blas/gemm.hpp"
 #include "oneapi/dal/backend/primitives/sparse_blas/gemv.hpp"
+#include "oneapi/dal/backend/primitives/sparse_blas/matmatd.hpp"

cpp/oneapi/dal/backend/primitives/sparse_blas/matmatd.hpp

@@ -0,0 +1,69 @@
+/*******************************************************************************
+* Copyright contributors to the oneDAL project
+*
+* Licensed under the Apache License, Version 2.0 (the "License");
+* you may not use this file except in compliance with the License.
+* You may obtain a copy of the License at
+*
+*     http://www.apache.org/licenses/LICENSE-2.0
+*
+* Unless required by applicable law or agreed to in writing, software
+* distributed under the License is distributed on an "AS IS" BASIS,
+* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+* See the License for the specific language governing permissions and
+* limitations under the License.
+*******************************************************************************/
+
+#pragma once
+
+#include "oneapi/dal/backend/primitives/ndarray.hpp"
+#include "oneapi/dal/backend/primitives/sparse_blas/misc.hpp"
+#include "oneapi/dal/backend/primitives/sparse_blas/handle.hpp"
+
+namespace oneapi::dal::backend::primitives {
+
+#ifdef ONEDAL_DATA_PARALLEL
+
+/// Computes a sparse matrix - sparse matrix product with a dense result:
+///         C = alpha * op(A) * op(B) + beta * C
+/// where `alpha` and `beta` are scalars, A, B - sparse matrices and C - dense matrix.
+/// op(A), op(B) are operators defining if the matrices A and B used as is in the computations
+/// or is being transposed.
+///
+/// op(A) is `m` x `k` matrix;
+/// op(B) is `k` x `p` matrix;
+/// C is `m` x `p` matrix.
+///
+/// @tparam Float   The type of elements in the matrix A and the vectors x and y.
+///                 The `Float` type should be at least `float` or `double`.
+/// @tparam co      Data layout in the matrix C.
+///                 The `co` value should be `ndorder::c` or `ndorder::f`.
+///
+/// @param queue        The SYCL* queue object.
+/// @param transpose_a  Defines if the sparse matrix A transposed or not.
+///                     If `transpose_a` == `transpose::notrans` then op(A) = A.
+///                     If `transpose_a` == `transpose::trans` then op(A) = transpose(A).
+/// @param transpose_b  Defines if the sparse matrix B transposed or not.
+///                     If `transpose_b` == `transpose::notrans` then op(B) =B.
+///                     If `transpose_b` == `transpose::trans` then op(B) = transpose(B).
+/// @param alpha        Specifies the scalar `alpha`.
+/// @param a            Handle to object containing sparse matrix A.
+/// @param b            Handle to object containing sparse matrix B.
+/// @param beta         Specifies the scalar `beta`.
+/// @param c            Dense output matrix that has `m` rows and `p` columns.
+/// @param dependencies Events indicating availability of the matrix A and the vectors x and y
+///                     for reading or writing.
+template <typename Float, ndorder co>
+sycl::event matmatd(sycl::queue &queue,
+                    transpose transpose_a,
+                    transpose transpose_b,
+                    const Float alpha,
+                    sparse_matrix_handle& a,
+                    sparse_matrix_handle& b,
+                    const Float beta,
+                    ndview<Float, 2, co>& c,
+                    const std::vector<sycl::event> &dependencies = {});
+
+#endif // ifdef ONEDAL_DATA_PARALLEL
+
+} // namespace oneapi::dal::backend::primitives

cpp/oneapi/dal/backend/primitives/sparse_blas/matmatd_dpc.cpp

@@ -0,0 +1,71 @@
+/*******************************************************************************
+* Copyright contributors to the oneDAL project
+*
+* Licensed under the Apache License, Version 2.0 (the "License");
+* you may not use this file except in compliance with the License.
+* You may obtain a copy of the License at
+*
+*     http://www.apache.org/licenses/LICENSE-2.0
+*
+* Unless required by applicable law or agreed to in writing, software
+* distributed under the License is distributed on an "AS IS" BASIS,
+* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+* See the License for the specific language governing permissions and
+* limitations under the License.
+*******************************************************************************/
+
+#include "oneapi/dal/backend/primitives/blas/misc.hpp"
+#include "oneapi/dal/backend/primitives/sparse_blas/matmatd.hpp"
+
+namespace oneapi::dal::backend::primitives {
+
+template <typename Float, ndorder co>
+sycl::event matmatd(sycl::queue &queue,
+                    transpose transpose_a,
+                    transpose transpose_b,
+                    const Float alpha,
+                    sparse_matrix_handle& a,
+                    sparse_matrix_handle& b,
+                    const Float beta,
+                    ndview<Float, 2, co>& c,
+                    const std::vector<sycl::event> &dependencies) {
+    ONEDAL_ASSERT(c.has_mutable_data());
+    mkl::transpose mkl_trans_a = transpose_to_mkl(transpose_a);
+    mkl::transpose mkl_trans_b = transpose_to_mkl(transpose_b);
+
+    if (co == ndorder::c) {
+        return mkl::sparse::matmatd(queue,
+                                    order_as_layout(co),
+                                    mkl_trans_a,
+                                    mkl_trans_b,
+                                    alpha,
+                                    dal::detail::get_impl(a).get(),
+                                    dal::detail::get_impl(b).get(),
+                                    beta,
+                                    c.get_mutable_data(),
+                                    c.get_dimension(0),
+                                    c.get_dimension(1),
+                                    c.get_leading_stride(),
+                                    dependencies);
+    }
+
+    return sycl::event();
+}
+
+#define INSTANTIATE(F, co)                                                        \
+    template ONEDAL_EXPORT sycl::event matmatd<F, co>(sycl::queue & queue,        \
+                                                      transpose transpose_a,      \
+                                                      transpose transpose_b,      \
+                                                      const F alpha,              \
+                                                      sparse_matrix_handle& a,    \
+                                                      sparse_matrix_handle& b,    \
+                                                      const F beta,               \
+                                                      ndview<F, 2, co>& c,        \
+                                                      const std::vector<sycl::event> &dependencies);
+
+INSTANTIATE(float, ndorder::c);
+INSTANTIATE(float, ndorder::f);
+INSTANTIATE(double, ndorder::c);
+INSTANTIATE(double, ndorder::f);
+
+} // namespace oneapi::dal::backend::primitives

cpp/oneapi/dal/backend/primitives/sparse_blas/test/fixture.hpp

@@ -85,13 +85,19 @@ class sparse_blas_test : public te::float_algo_fixture<std::tuple_element_t<0, P
     /// | 1  0  1  0  ... |
     /// | 0  1  0  1  ... |
     /// | ............... |
-    auto A(sparse_matrix_handle& a) {
+    auto A(sparse_matrix_handle& h,
+           const transpose& trans,
+           std::int64_t m,
+           std::int64_t k,
+           dal::array<float_t>& data_ary_,
+           dal::array<std::int64_t>& column_indices_ary_,
+           dal::array<std::int64_t>& row_offsets_ary_) {
         check_if_initialized();
         auto& q = this->get_queue();
 
         /// Revert dimensions in the transposed case
-        const std::int64_t local_m = trans_a == transpose::nontrans ? m_ : k_;
-        const std::int64_t local_k = trans_a == transpose::nontrans ? k_ : m_;
+        const std::int64_t local_m = trans == transpose::nontrans ? m : k;
+        const std::int64_t local_k = trans == transpose::nontrans ? k : m;
 
         const std::int64_t value_count_in_odd_rows = (local_k + 1) / 2;
         const std::int64_t value_count_in_even_rows = local_k / 2;
@@ -157,7 +163,7 @@ class sparse_blas_test : public te::float_algo_fixture<std::tuple_element_t<0, P
         });
 
         return set_csr_data(q,
-                            a,
+                            h,
                             local_m,
                             local_k,
                             indexing,
@@ -185,7 +191,7 @@ class sparse_blas_test : public te::float_algo_fixture<std::tuple_element_t<0, P
 
     void test_gemm() {
         sparse_matrix_handle a(this->get_queue());
-        auto a_e = A(a);
+        auto a_e = A(a, trans_a, m_, k_, a_data_ary_, a_column_indices_ary_, a_row_offsets_ary_);
         auto [b, b_e] = B();
         auto c = C();
 
@@ -194,9 +200,22 @@ class sparse_blas_test : public te::float_algo_fixture<std::tuple_element_t<0, P
         check_matmul(c);
     }
 
+    void test_matmatd() {
+        sparse_matrix_handle a(this->get_queue());
+        sparse_matrix_handle b(this->get_queue());
+        constexpr transpose trans_b  = (bo == ndorder::c ? transpose::nontrans : transpose::trans);
+        auto a_e = A(a, trans_a, m_, k_, a_data_ary_, a_column_indices_ary_, a_row_offsets_ary_);
+        auto b_e = A(b, trans_b, k_, p_, b_data_ary_, b_column_indices_ary_, b_row_offsets_ary_);
+        auto c = C();
+
+        matmatd(this->get_queue(), trans_a, trans_b, float_t(1.0), a, b, float_t(0.0), c, { a_e, b_e }).wait_and_throw();
+
+        check_matmul(c);
+    }
+
     void test_gemv() {
         sparse_matrix_handle a(this->get_queue());
-        auto a_e = A(a);
+        auto a_e = A(a, trans_a, m_, k_, a_data_ary_, a_column_indices_ary_, a_row_offsets_ary_);
         auto [x, x_e] = ndarray<float_t, 1>::ones(this->get_queue(), k_);
         auto y = ndarray<float_t, 1>::empty(this->get_queue(), m_);
 
@@ -269,9 +288,14 @@ class sparse_blas_test : public te::float_algo_fixture<std::tuple_element_t<0, P
     std::int64_t p_;
 
     /// Sparse matrix A
-    dal::array<float_t> data_ary_;
-    dal::array<std::int64_t> column_indices_ary_;
-    dal::array<std::int64_t> row_offsets_ary_;
+    dal::array<float_t> a_data_ary_;
+    dal::array<std::int64_t> a_column_indices_ary_;
+    dal::array<std::int64_t> a_row_offsets_ary_;
+
+    /// Sparse matrix B
+    dal::array<float_t> b_data_ary_;
+    dal::array<std::int64_t> b_column_indices_ary_;
+    dal::array<std::int64_t> b_row_offsets_ary_;
 
     sycl::usm::alloc alloc_;
 };

cpp/oneapi/dal/backend/primitives/sparse_blas/test/matmatd_dpc.cpp

@@ -0,0 +1,41 @@
+/*******************************************************************************
+* Copyright contributors to the oneDAL project
+*
+* Licensed under the Apache License, Version 2.0 (the "License");
+* you may not use this file except in compliance with the License.
+* You may obtain a copy of the License at
+*
+*     http://www.apache.org/licenses/LICENSE-2.0
+*
+* Unless required by applicable law or agreed to in writing, software
+* distributed under the License is distributed on an "AS IS" BASIS,
+* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+* See the License for the specific language governing permissions and
+* limitations under the License.
+*******************************************************************************/
+
+#include "oneapi/dal/backend/primitives/sparse_blas/test/fixture.hpp"
+
+namespace oneapi::dal::backend::primitives::test {
+
+namespace te = dal::test::engine;
+
+using matmatd_types = COMBINE_TYPES((float, double),
+                                    (transpose_nontrans, transpose_trans),
+                                    (c_order, f_order),
+                                    (c_order), /// not used in matmatd
+                                    (indexing_zero_based, indexing_one_based));
+
+TEMPLATE_LIST_TEST_M(sparse_blas_test, "ones matrix sparse CSR matmatd", "[csr][matmatd]", matmatd_types) {
+    // DPC++ Sparse GEMV from micro MKL libs is not supported on CPU
+    SKIP_IF(this->get_policy().is_cpu());
+
+    // Test takes too long time if HW emulates float64
+    // Temporary workaround: skip tests on architectures that do not support native float64
+    SKIP_IF(!this->get_policy().has_native_float64());
+
+    this->generate_dimensions();
+    this->test_matmatd();
+}
+
+} // namespace oneapi::dal::backend::primitives::test