EIG combine existing MKL and non-existing SYCL to one kernel

Author: densmirn

dpnp/backend/backend_iface.hpp

@@ -189,7 +189,7 @@ INP_DLLEXPORT void mkl_blas_dot_c(void* array1, void* array2, void* result1, siz
 
 /**
  * @ingroup BACKEND_API
- * @brief MKL implementation of eig function
+ * @brief Custom implementation of eig function
  *
  * @param [in]  array1  Input array.
  *
@@ -199,7 +199,7 @@ INP_DLLEXPORT void mkl_blas_dot_c(void* array1, void* array2, void* result1, siz
  *
  */
 template <typename _DataType>
-INP_DLLEXPORT void mkl_lapack_syevd_c(void* array1, void* result1, size_t size);
+INP_DLLEXPORT void custom_lapack_syevd_c(void* array1, void* result1, size_t size);
 
 /**
  * @ingroup BACKEND_API

dpnp/backend/backend_iface_fptr.cpp

@@ -315,10 +315,10 @@ static func_map_t func_map_init()
     fmap[DPNPFuncName::DPNP_FN_DOT][eft_FLT][eft_FLT] = {eft_FLT, (void*)mkl_blas_dot_c<float>};
     fmap[DPNPFuncName::DPNP_FN_DOT][eft_DBL][eft_DBL] = {eft_DBL, (void*)mkl_blas_dot_c<double>};
 
-    fmap[DPNPFuncName::DPNP_FN_EIG][eft_INT][eft_INT] = {eft_DBL, (void*)mkl_lapack_syevd_c<double>};
-    fmap[DPNPFuncName::DPNP_FN_EIG][eft_LNG][eft_LNG] = {eft_DBL, (void*)mkl_lapack_syevd_c<double>};
-    fmap[DPNPFuncName::DPNP_FN_EIG][eft_FLT][eft_FLT] = {eft_FLT, (void*)mkl_lapack_syevd_c<float>};
-    fmap[DPNPFuncName::DPNP_FN_EIG][eft_DBL][eft_DBL] = {eft_DBL, (void*)mkl_lapack_syevd_c<double>};
+    fmap[DPNPFuncName::DPNP_FN_EIG][eft_INT][eft_INT] = {eft_DBL, (void*)custom_lapack_syevd_c<double>};
+    fmap[DPNPFuncName::DPNP_FN_EIG][eft_LNG][eft_LNG] = {eft_DBL, (void*)custom_lapack_syevd_c<double>};
+    fmap[DPNPFuncName::DPNP_FN_EIG][eft_FLT][eft_FLT] = {eft_FLT, (void*)custom_lapack_syevd_c<float>};
+    fmap[DPNPFuncName::DPNP_FN_EIG][eft_DBL][eft_DBL] = {eft_DBL, (void*)custom_lapack_syevd_c<double>};
 
     fmap[DPNPFuncName::DPNP_FN_EXP][eft_INT][eft_INT] = {eft_DBL, (void*)custom_elemwise_exp_c<int, double>};
     fmap[DPNPFuncName::DPNP_FN_EXP][eft_LNG][eft_LNG] = {eft_DBL, (void*)custom_elemwise_exp_c<long, double>};

dpnp/backend/custom_kernels.cpp

@@ -26,6 +26,7 @@
 #include <cmath>
 #include <iostream>
 #include <mkl_blas_sycl.hpp>
+#include <mkl_lapack_sycl.hpp>
 #include <type_traits>
 
 #include <backend_iface.hpp>
@@ -34,6 +35,7 @@
 #include "queue_sycl.hpp"
 
 namespace mkl_blas = oneapi::mkl::blas;
+namespace mkl_lapack = oneapi::mkl::lapack;
 
 template <typename _KernelNameSpecialization>
 class custom_blas_gemm_c_kernel;
@@ -161,6 +163,60 @@ void custom_blas_dot_c(void* array1_in, void* array2_in, void* result1, size_t s
 template void custom_blas_dot_c<long>(void* array1_in, void* array2_in, void* result1, size_t size);
 template void custom_blas_dot_c<int>(void* array1_in, void* array2_in, void* result1, size_t size);
 
+template <typename _DataType>
+void custom_lapack_syevd_c(void* array_in, void* result1, size_t size)
+{
+    if (!size)
+    {
+        return;
+    }
+
+    _DataType* array = reinterpret_cast<_DataType*>(array_in);
+    _DataType* result = reinterpret_cast<_DataType*>(result1);
+
+    if constexpr (std::is_same<_DataType, double>::value || std::is_same<_DataType, float>::value)
+    {
+        cl::sycl::event event;
+
+        _DataType* syevd_array = reinterpret_cast<_DataType*>(dpnp_memory_alloc_c(size * size * sizeof(_DataType)));
+        dpnp_memory_memcpy_c(syevd_array, array, size * size * sizeof(_DataType));
+
+        const std::int64_t lda = std::max<size_t>(1UL, size);
+
+        const std::int64_t scratchpad_size = mkl_lapack::syevd_scratchpad_size<_DataType>(
+            DPNP_QUEUE, oneapi::mkl::job::vec, oneapi::mkl::uplo::upper, size, lda);
+
+        _DataType* scratchpad = reinterpret_cast<_DataType*>(dpnp_memory_alloc_c(scratchpad_size * sizeof(_DataType)));
+
+        event = mkl_lapack::syevd(DPNP_QUEUE,               // queue
+                                  oneapi::mkl::job::vec,    // jobz
+                                  oneapi::mkl::uplo::upper, // uplo
+                                  size,                     // The order of the matrix A (0≤n)
+                                  syevd_array,              // will be overwritten with eigenvectors
+                                  lda,
+                                  result,
+                                  scratchpad,
+                                  scratchpad_size);
+        event.wait();
+
+        dpnp_memory_free_c(scratchpad);
+
+        custom_elemwise_transpose_c<_DataType>(
+            syevd_array, {(long)size, (long)size}, {(long)size, (long)size}, {1, 0}, array, size * size);
+
+        dpnp_memory_free_c(syevd_array);
+    }
+    else
+    {
+        // TODO: implement SYCL kernel for int/long input
+    }
+}
+
+template void custom_lapack_syevd_c<int>(void* array1_in, void* result1, size_t size);
+template void custom_lapack_syevd_c<long>(void* array1_in, void* result1, size_t size);
+template void custom_lapack_syevd_c<float>(void* array1_in, void* result1, size_t size);
+template void custom_lapack_syevd_c<double>(void* array1_in, void* result1, size_t size);
+
 #if 0 // Example for OpenCL kernel
 #include <map>
 #include <typeindex>

dpnp/backend/mkl_wrap_lapack.cpp

@@ -28,7 +28,7 @@ int main(int, char**)
         array[size * i + i] = i + 1;
     }
 
-    mkl_lapack_syevd_c<double>(array, result, size);
+    custom_lapack_syevd_c<double>(array, result, size);
 
     std::cout << "eigen values" << std::endl;
     for (size_t i = 0; i < size; ++i)

examples/example7.cpp

@@ -273,7 +273,6 @@
                 "dpnp/backend/custom_kernels_statistics.cpp",
                 "dpnp/backend/memory_sycl.cpp",
                 "dpnp/backend/mkl_wrap_blas1.cpp",
-                "dpnp/backend/mkl_wrap_lapack.cpp",
                 "dpnp/backend/mkl_wrap_rng.cpp",
                 "dpnp/backend/queue_sycl.cpp"
             ],