linalg.eig and beta09 related changes

dpnp/backend/mkl_wrap_blas1.cpp

@@ -30,6 +30,8 @@
 #include <backend/backend_iface.hpp>
 #include "queue_sycl.hpp"
 
+namespace mkl_blas = oneapi::mkl::blas;
+
 template <typename _DataType>
 void mkl_blas_dot_c(void* array1_in, void* array2_in, void* result1, size_t size)
 {
@@ -41,13 +43,13 @@ void mkl_blas_dot_c(void* array1_in, void* array2_in, void* result1, size_t size
 
     try
     {
-        status = mkl::blas::dot(DPNP_QUEUE,
-                                size,
-                                array_1,
-                                1, // array_1 stride
-                                array_2,
-                                1, // array_2 stride
-                                result);
+        status = mkl_blas::dot(DPNP_QUEUE,
+                               size,
+                               array_1,
+                               1, // array_1 stride
+                               array_2,
+                               1, // array_2 stride
+                               result);
     }
     catch (cl::sycl::exception const& e)
     {

dpnp/backend/mkl_wrap_blas3.cpp

@@ -30,6 +30,8 @@
 #include <backend/backend_iface.hpp>
 #include "queue_sycl.hpp"
 
+namespace mkl_blas = oneapi::mkl::blas;
+
 template <typename _DataType>
 void dpnp_blas_gemm_c(void* array1_in, void* array2_in, void* result1, size_t size_m, size_t size_n, size_t size_k)
 {
@@ -89,20 +91,20 @@ void dpnp_blas_gemm_c(void* array1_in, void* array2_in, void* result1, size_t si
 #endif
     try
     {
-        status = mkl::blas::gemm(DPNP_QUEUE,
-                                 mkl::transpose::nontrans,
-                                 mkl::transpose::nontrans,
-                                 size_n,
-                                 size_m,
-                                 size_k,
-                                 _DataType(1),
-                                 array_2,
-                                 ldb,
-                                 array_1,
-                                 lda,
-                                 _DataType(0),
-                                 result,
-                                 ldc);
+        status = mkl_blas::gemm(DPNP_QUEUE,
+                                oneapi::mkl::transpose::nontrans,
+                                oneapi::mkl::transpose::nontrans,
+                                size_n,
+                                size_m,
+                                size_k,
+                                _DataType(1),
+                                array_2,
+                                ldb,
+                                array_1,
+                                lda,
+                                _DataType(0),
+                                result,
+                                ldc);
     }
     catch (cl::sycl::exception const& e)
     {

dpnp/backend/mkl_wrap_lapack.cpp

@@ -30,6 +30,8 @@
 #include <backend/backend_iface.hpp>
 #include "queue_sycl.hpp"
 
+namespace mkl_lapack = oneapi::mkl::lapack;
+
 template <typename _DataType>
 void mkl_lapack_syevd_c(void* array_in, void* result1, size_t size)
 {
@@ -38,60 +40,33 @@ void mkl_lapack_syevd_c(void* array_in, void* result1, size_t size)
     _DataType* array = reinterpret_cast<_DataType*>(array_in);
     _DataType* result = reinterpret_cast<_DataType*>(result1);
 
-    const std::int64_t lda = std::max<size_t>(1UL, size);
+    _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));
 
-    auto queue = DPNP_QUEUE;
+    const std::int64_t lda = std::max<size_t>(1UL, size);
 
-    const std::int64_t scratchpad_size =
-        mkl::lapack::syevd_scratchpad_size<_DataType>(queue, mkl::job::vec, mkl::uplo::upper, size, lda);
-    // const std::int64_t scratchpad_size = 7; // size = 2 & novec & lower
-    // const std::int64_t scratchpad_size = 34; // size = 2 & vec & upper
+    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)));
 
-#if 1
-    printf("mkl_lapack_syevd_c array ptr type %d \n", cl::sycl::get_pointer_type(array, queue.get_context()));
-    printf("mkl_lapack_syevd_c result ptr type %d \n", cl::sycl::get_pointer_type(result, queue.get_context()));
-    printf("mkl_lapack_syevd_c scratchpad ptr type %d \n", cl::sycl::get_pointer_type(scratchpad, queue.get_context()));
-    printf("mkl_lapack_syevd_c array[0] %g \n", array[0]);
-    printf("mkl_lapack_syevd_c result[0] %g \n", result[0]);
-    std::cout << "mkl_lapack_syevd_c size = " << size << std::endl;
-    std::cout << "mkl_lapack_syevd_c lda = " << lda << std::endl;
-    std::cout << "mkl_lapack_syevd_c scratchpad_size = " << scratchpad_size << std::endl;
-#endif
-
-    try
-    {
-        queue.wait_and_throw();
-        status = mkl::lapack::syevd(queue,            // queue
-                                    mkl::job::vec,    // jobz
-                                    mkl::uplo::upper, // uplo
-                                    size,             // The order of the matrix A (0≤n)
-                                    array,            // will be overwritten with eigenvectors
-                                    lda,
-                                    result,
-                                    scratchpad,
-                                    scratchpad_size);
-    }
-    catch (mkl::lapack::exception const& e)
-    {
-        // Handle LAPACK related exceptions happened during asynchronous call
-        std::cout << "Unexpected exception caught during asynchronous LAPACK operation:\n"
-                  << e.reason() << "\ninfo: " << e.info() << std::endl;
-    }
-    catch (cl::sycl::exception const& e)
-    {
-        std::cerr << "Caught synchronous SYCL exception during mkl_lapack_syevd_c():\n"
-                  << e.what() << "\nOpenCL status: " << e.get_cl_code() << std::endl;
-    }
-
+    status = 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);
     status.wait();
 
-#if 1
-    std::cout << "mkl_lapack_syevd_c res = " << result[0] << std::endl;
-#endif
-
     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);
 }
 
 template void mkl_lapack_syevd_c<double>(void* array1_in, void* result1, size_t size);

dpnp/backend/mkl_wrap_rng.cpp

@@ -33,6 +33,8 @@
 #include "backend_utils.hpp"
 #include "queue_sycl.hpp"
 
+namespace mkl_rng = oneapi::mkl::rng;
+
 template <typename _DataType>
 void mkl_rng_gaussian(void* result, size_t size)
 {
@@ -42,16 +44,16 @@ void mkl_rng_gaussian(void* result, size_t size)
     // choose engine as is in numpy
     // seed number
     size_t seed = std::time(nullptr);
-    mkl::rng::philox4x32x10 engine(DPNP_QUEUE, seed);
+    mkl_rng::philox4x32x10 engine(DPNP_QUEUE, seed);
 
     const _DataType mean = _DataType(0.0);
     const _DataType stddev = _DataType(1.0);
 
-    mkl::rng::gaussian<_DataType> distribution(mean, stddev);
+    mkl_rng::gaussian<_DataType> distribution(mean, stddev);
     try
     {
         // perform generation
-        mkl::rng::generate(distribution, engine, size, result1);
+        mkl_rng::generate(distribution, engine, size, result1);
         DPNP_QUEUE.wait_and_throw();
     }
     catch (cl::sycl::exception const& e)
@@ -70,16 +72,16 @@ void mkl_rng_uniform(void* result, size_t size)
     // choose engine as is in numpy
     // seed number
     size_t seed = std::time(nullptr);
-    mkl::rng::philox4x32x10 engine(DPNP_QUEUE, seed);
+    mkl_rng::philox4x32x10 engine(DPNP_QUEUE, seed);
 
     const _DataType a = (_DataType(0.0));
     const _DataType b = (_DataType(1.0));
 
-    mkl::rng::uniform<_DataType> distribution(a, b);
+    mkl_rng::uniform<_DataType> distribution(a, b);
     try
     {
         // perform generation
-        mkl::rng::generate(distribution, engine, size, result1);
+        mkl_rng::generate(distribution, engine, size, result1);
         DPNP_QUEUE.wait_and_throw();
     }
     catch (cl::sycl::exception const& e)
@@ -98,18 +100,18 @@ void mkl_rng_uniform_mt19937(void* result, long low, long high, size_t size)
     // choose engine as is in numpy
     // seed number
     size_t seed = std::time(nullptr);
-    mkl::rng::mt19937 engine(DPNP_QUEUE, seed);
+    mkl_rng::mt19937 engine(DPNP_QUEUE, seed);
 
     // set left bound of distribution
     const _DataType a = (_DataType(low));
     // set right bound of distribution
     const _DataType b = (_DataType(high));
 
-    mkl::rng::uniform<_DataType> distribution(a, b);
+    mkl_rng::uniform<_DataType> distribution(a, b);
     try
     {
         // perform generation
-        mkl::rng::generate(distribution, engine, size, result1);
+        mkl_rng::generate(distribution, engine, size, result1);
         DPNP_QUEUE.wait_and_throw();
     }
     catch (cl::sycl::exception const& e)

dpnp/linalg/linalg.pyx

@@ -45,23 +45,26 @@ __all__ = [
 ]
 
 
-cpdef dparray dpnp_eig(dparray in_array1):
-    cdef vector[Py_ssize_t] shape1 = in_array1.shape
+cpdef tuple dpnp_eig(dparray in_array1):
+    cdef dparray_shape_type shape1 = in_array1.shape
 
     call_type = in_array1.dtype
+    res_type = call_type
+    if res_type != numpy.float32:
+        res_type = numpy.float64
 
     cdef size_t size1 = 0
     if not shape1.empty():
         size1 = shape1.front()
 
-    cdef dparray res_val = dparray((size1,), dtype=call_type)
+    cdef dparray res_val = dparray((size1,), dtype=res_type)
 
     # this array is used as input for MKL and will be overwritten with eigen vectors
-    cdef dparray res_vec = dparray(shape1, dtype=call_type)
+    cdef dparray res_vec = dparray(shape1, dtype=res_type)
     for i in range(in_array1.size):
         res_vec[i] = in_array1[i]
 
-    if call_type == numpy.float64:
+    if call_type in [numpy.float64, numpy.int32, numpy.int64]:
         mkl_lapack_syevd_c[double](res_vec.get_data(), res_val.get_data(), size1)
     elif call_type == numpy.float32:
         mkl_lapack_syevd_c[float](res_vec.get_data(), res_val.get_data(), size1)

tests/test_linalg.py

@@ -0,0 +1,50 @@
+import pytest
+
+import dpnp as inp
+
+import numpy
+
+
+def vvsort(val, vec, size):
+    for i in range(size):
+        imax = i
+        for j in range(i + 1, size):
+            if numpy.abs(val[imax]) < numpy.abs(val[j]):
+                imax = j
+
+        temp = val[i]
+        val[i] = val[imax]
+        val[imax] = temp
+
+        for k in range(size):
+            temp = vec[k, i]
+            vec[k, i] = vec[k, imax]
+            vec[k, imax] = temp
+
+
+@pytest.mark.parametrize("type",
+                         [numpy.float64, numpy.float32, numpy.int64, numpy.int32],
+                         ids=['float64', 'float32', 'int64', 'int32'])
+@pytest.mark.parametrize("size",
+                         [2, 4, 8, 16, 300])
+def test_eig_arange(type, size):
+    a = numpy.arange(size * size, dtype=type).reshape((size, size))
+    symm = numpy.tril(a) + numpy.tril(a, -1).T + numpy.diag(numpy.full((size,), size * size, dtype=type))
+    isymm = inp.array(symm)
+
+    dpnp_val, dpnp_vec = inp.linalg.eig(isymm)
+    np_val, np_vec = numpy.linalg.eig(symm)
+
+    # DPNP sort val/vec by abs value
+    vvsort(dpnp_val, dpnp_vec, size)
+
+    # NP sort val/vec by abs value
+    vvsort(np_val, np_vec, size)
+
+    # NP change sign of vectors
+    for i in range(np_vec.shape[1]):
+        if np_vec[0, i] * dpnp_vec[0, i] < 0:
+            np_vec[:, i] = -np_vec[:, i]
+
+    numpy.testing.assert_allclose(dpnp_val, np_val, rtol=1e-05, atol=1e-05)
+    numpy.testing.assert_allclose(dpnp_vec, np_vec, rtol=1e-05, atol=1e-05)