diff --git a/dpnp/backend/mkl_wrap_blas1.cpp b/dpnp/backend/mkl_wrap_blas1.cpp index 230d7e8253fd..635f1ae73651 100644 --- a/dpnp/backend/mkl_wrap_blas1.cpp +++ b/dpnp/backend/mkl_wrap_blas1.cpp @@ -30,6 +30,8 @@ #include #include "queue_sycl.hpp" +namespace mkl_blas = oneapi::mkl::blas; + template 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) { diff --git a/dpnp/backend/mkl_wrap_blas3.cpp b/dpnp/backend/mkl_wrap_blas3.cpp index 6e21883a6131..5bc69abe1530 100644 --- a/dpnp/backend/mkl_wrap_blas3.cpp +++ b/dpnp/backend/mkl_wrap_blas3.cpp @@ -30,6 +30,8 @@ #include #include "queue_sycl.hpp" +namespace mkl_blas = oneapi::mkl::blas; + template 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) { diff --git a/dpnp/backend/mkl_wrap_lapack.cpp b/dpnp/backend/mkl_wrap_lapack.cpp index 64c24f72b5cb..e88c36305ea3 100644 --- a/dpnp/backend/mkl_wrap_lapack.cpp +++ b/dpnp/backend/mkl_wrap_lapack.cpp @@ -30,6 +30,8 @@ #include #include "queue_sycl.hpp" +namespace mkl_lapack = oneapi::mkl::lapack; + template 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(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(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(void* array1_in, void* result1, size_t size); diff --git a/dpnp/backend/mkl_wrap_rng.cpp b/dpnp/backend/mkl_wrap_rng.cpp index ce89b38741a0..ffdb8fb3009e 100644 --- a/dpnp/backend/mkl_wrap_rng.cpp +++ b/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 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) diff --git a/dpnp/linalg/linalg.pyx b/dpnp/linalg/linalg.pyx index f5fa11c19be7..3d4647263cef 100644 --- a/dpnp/linalg/linalg.pyx +++ b/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) diff --git a/tests/test_linalg.py b/tests/test_linalg.py new file mode 100644 index 000000000000..13b188699d43 --- /dev/null +++ b/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)