Merge DPCPP SpGEMM, SpGEAM, Transpose, Sort

This PR adds CSR sparse BLAS kernels (SpGEMM, SpGEAM, Transpose, Sort) with
naive implementations to provide the necessary functionality.

Related PR: #799

dpcpp/components/atomic.dp.hpp

@@ -0,0 +1,211 @@
+/*******************************<GINKGO LICENSE>******************************
+Copyright (c) 2017-2021, the Ginkgo authors
+All rights reserved.
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions
+are met:
+
+1. Redistributions of source code must retain the above copyright
+notice, this list of conditions and the following disclaimer.
+
+2. Redistributions in binary form must reproduce the above copyright
+notice, this list of conditions and the following disclaimer in the
+documentation and/or other materials provided with the distribution.
+
+3. Neither the name of the copyright holder nor the names of its
+contributors may be used to endorse or promote products derived from
+this software without specific prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS
+IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
+TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
+PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
+HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
+SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
+LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
+DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
+THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+******************************<GINKGO LICENSE>*******************************/
+
+#ifndef GKO_DPCPP_COMPONENTS_ATOMIC_DP_HPP_
+#define GKO_DPCPP_COMPONENTS_ATOMIC_DP_HPP_
+
+
+#include <type_traits>
+
+
+#include <CL/sycl.hpp>
+
+
+#include "dpcpp/base/dpct.hpp"
+
+
+namespace gko {
+namespace kernels {
+namespace dpcpp {
+namespace atomic {
+
+
+constexpr auto local_space = cl::sycl::access::address_space::local_space;
+constexpr auto global_space = cl::sycl::access::address_space::global_space;
+
+
+}  // namespace atomic
+
+namespace {
+
+
+template <cl::sycl::access::address_space addressSpace = atomic::global_space,
+          typename T>
+T atomic_compare_exchange_strong(
+    cl::sycl::multi_ptr<T, addressSpace> addr, T expected, T desired,
+    cl::sycl::memory_order success = cl::sycl::memory_order::relaxed,
+    cl::sycl::memory_order fail = cl::sycl::memory_order::relaxed)
+{
+    cl::sycl::atomic<T, addressSpace> obj(addr);
+    obj.compare_exchange_strong(expected, desired, success, fail);
+    return expected;
+}
+
+template <cl::sycl::access::address_space addressSpace = atomic::global_space,
+          typename T>
+T atomic_compare_exchange_strong(
+    T *addr, T expected, T desired,
+    cl::sycl::memory_order success = cl::sycl::memory_order::relaxed,
+    cl::sycl::memory_order fail = cl::sycl::memory_order::relaxed)
+{
+    return atomic_compare_exchange_strong(
+        cl::sycl::multi_ptr<T, addressSpace>(addr), expected, desired, success,
+        fail);
+}
+
+
+template <cl::sycl::access::address_space addressSpace = atomic::global_space,
+          typename T>
+inline T atomic_fetch_add(
+    T *addr, T operand,
+    cl::sycl::memory_order memoryOrder = cl::sycl::memory_order::relaxed)
+{
+    cl::sycl::atomic<T, addressSpace> obj(
+        (cl::sycl::multi_ptr<T, addressSpace>(addr)));
+    return cl::sycl::atomic_fetch_add(obj, operand, memoryOrder);
+}
+
+
+}  // namespace
+
+
+namespace detail {
+
+
+template <cl::sycl::access::address_space addressSpace, typename ValueType,
+          typename = void>
+struct atomic_helper {
+    __dpct_inline__ static ValueType atomic_add(ValueType *, ValueType)
+    {
+        static_assert(sizeof(ValueType) == 0,
+                      "This default function is not implemented, only the "
+                      "specializations are.");
+        // TODO: add proper implementation of generic atomic add
+    }
+};
+
+
+template <typename ResultType, typename ValueType>
+__dpct_inline__ ResultType reinterpret(ValueType val)
+{
+    static_assert(sizeof(ValueType) == sizeof(ResultType),
+                  "The type to reinterpret to must be of the same size as the "
+                  "original type.");
+    return reinterpret_cast<ResultType &>(val);
+}
+
+
+#define GKO_BIND_ATOMIC_HELPER_STRUCTURE(CONVERTER_TYPE)                   \
+    template <cl::sycl::access::address_space addressSpace,                \
+              typename ValueType>                                          \
+    struct atomic_helper<                                                  \
+        addressSpace, ValueType,                                           \
+        std::enable_if_t<(sizeof(ValueType) == sizeof(CONVERTER_TYPE))>> { \
+        __dpct_inline__ static ValueType atomic_add(                       \
+            ValueType *__restrict__ addr, ValueType val)                   \
+        {                                                                  \
+            CONVERTER_TYPE *address_as_converter =                         \
+                reinterpret_cast<CONVERTER_TYPE *>(addr);                  \
+            CONVERTER_TYPE old = *address_as_converter;                    \
+            CONVERTER_TYPE assumed;                                        \
+            do {                                                           \
+                assumed = old;                                             \
+                old = atomic_compare_exchange_strong<addressSpace>(        \
+                    address_as_converter, assumed,                         \
+                    reinterpret<CONVERTER_TYPE>(                           \
+                        val + reinterpret<ValueType>(assumed)));           \
+            } while (assumed != old);                                      \
+            return reinterpret<ValueType>(old);                            \
+        }                                                                  \
+    };
+
+// Support 64-bit ATOMIC_ADD
+GKO_BIND_ATOMIC_HELPER_STRUCTURE(unsigned long long int);
+// Support 32-bit ATOMIC_ADD
+GKO_BIND_ATOMIC_HELPER_STRUCTURE(unsigned int);
+
+
+#undef GKO_BIND_ATOMIC_HELPER_STRUCTURE
+
+#define GKO_BIND_ATOMIC_HELPER_VALUETYPE(ValueType)                         \
+    template <cl::sycl::access::address_space addressSpace>                 \
+    struct atomic_helper<addressSpace, ValueType, std::enable_if_t<true>> { \
+        __dpct_inline__ static ValueType atomic_add(                        \
+            ValueType *__restrict__ addr, ValueType val)                    \
+        {                                                                   \
+            return atomic_fetch_add<addressSpace>(addr, val);               \
+        }                                                                   \
+    };
+
+GKO_BIND_ATOMIC_HELPER_VALUETYPE(int);
+GKO_BIND_ATOMIC_HELPER_VALUETYPE(unsigned int);
+GKO_BIND_ATOMIC_HELPER_VALUETYPE(unsigned long long int);
+
+#undef GKO_BIND_ATOMIC_HELPER_VALUETYPE
+
+
+template <cl::sycl::access::address_space addressSpace, typename ValueType>
+struct atomic_helper<
+    addressSpace, ValueType,
+    std::enable_if_t<is_complex<ValueType>() && sizeof(ValueType) >= 16>> {
+    __dpct_inline__ static ValueType atomic_add(ValueType *__restrict__ addr,
+                                                ValueType val)
+    {
+        using real_type = remove_complex<ValueType>;
+        real_type *real_addr = reinterpret_cast<real_type *>(addr);
+        // Separate to real part and imag part
+        auto real = atomic_helper<addressSpace, real_type>::atomic_add(
+            &real_addr[0], val.real());
+        auto imag = atomic_helper<addressSpace, real_type>::atomic_add(
+            &real_addr[1], val.imag());
+        return {real, imag};
+    }
+};
+
+
+}  // namespace detail
+
+
+template <cl::sycl::access::address_space addressSpace = atomic::global_space,
+          typename T>
+__dpct_inline__ T atomic_add(T *__restrict__ addr, T val)
+{
+    return detail::atomic_helper<addressSpace, T>::atomic_add(addr, val);
+}
+
+
+}  // namespace dpcpp
+}  // namespace kernels
+}  // namespace gko
+
+
+#endif  // GKO_DPCPP_COMPONENTS_ATOMIC_DP_HPP_

dpcpp/components/prefix_sum.dp.cpp

@@ -47,7 +47,18 @@ namespace components {
 
 template <typename IndexType>
 void prefix_sum(std::shared_ptr<const DefaultExecutor> exec, IndexType *counts,
-                size_type num_entries) GKO_NOT_IMPLEMENTED;
+                size_type num_entries)
+{
+    // TODO actually implement parallel prefix sum
+    exec->get_queue()->submit([&](sycl::handler &cgh) {
+        cgh.parallel_for(sycl::range<1>{1}, [=](sycl::id<1> idx) {
+            IndexType sum{};
+            for (size_type i = 0; i < num_entries; i++) {
+                sum += std::exchange(counts[i], sum);
+            }
+        });
+    });
+}
 
 GKO_INSTANTIATE_FOR_EACH_INDEX_TYPE(GKO_DECLARE_PREFIX_SUM_KERNEL);
 // instantiate for size_type as well, as this is used in the Sellp format