From 57b92e219c901a1873ca02a905439b6ff0d8757c Mon Sep 17 00:00:00 2001 From: Anton Volkov Date: Wed, 19 Apr 2023 07:44:06 -0500 Subject: [PATCH 1/3] Add dpnp.linalg.eigh() function --- .gitignore | 1 + dpnp/CMakeLists.txt | 1 + dpnp/backend/CMakeLists.txt | 2 +- dpnp/backend/extensions/lapack/CMakeLists.txt | 75 ++++++ dpnp/backend/extensions/lapack/heevd.cpp | 210 +++++++++++++++++ dpnp/backend/extensions/lapack/heevd.hpp | 51 +++++ dpnp/backend/extensions/lapack/lapack_py.cpp | 57 +++++ dpnp/backend/extensions/lapack/syevd.cpp | 214 ++++++++++++++++++ dpnp/backend/extensions/lapack/syevd.hpp | 51 +++++ dpnp/dpnp_iface.py | 45 +++- dpnp/linalg/dpnp_iface_linalg.py | 67 ++++++ dpnp/linalg/dpnp_utils_linalg.py | 134 +++++++++++ tests/skipped_tests_gpu.tbl | 12 + tests/test_linalg.py | 47 +++- tests/test_sycl_queue.py | 32 +++ .../cupy/creation_tests/test_ranges.py | 2 +- .../cupy/linalg_tests/test_eigenvalue.py | 199 ++++++++++++++++ tests/third_party/cupy/testing/helper.py | 28 ++- 18 files changed, 1216 insertions(+), 12 deletions(-) create mode 100644 dpnp/backend/extensions/lapack/CMakeLists.txt create mode 100644 dpnp/backend/extensions/lapack/heevd.cpp create mode 100644 dpnp/backend/extensions/lapack/heevd.hpp create mode 100644 dpnp/backend/extensions/lapack/lapack_py.cpp create mode 100644 dpnp/backend/extensions/lapack/syevd.cpp create mode 100644 dpnp/backend/extensions/lapack/syevd.hpp create mode 100644 dpnp/linalg/dpnp_utils_linalg.py create mode 100644 tests/third_party/cupy/linalg_tests/test_eigenvalue.py diff --git a/.gitignore b/.gitignore index 8beb38f1efd6..ea9f2cba333d 100644 --- a/.gitignore +++ b/.gitignore @@ -1,5 +1,6 @@ # CMake build and local install directory _skbuild +build build_cython dpnp.egg-info diff --git a/dpnp/CMakeLists.txt b/dpnp/CMakeLists.txt index 8309b0422e4e..54be4eb23b9a 100644 --- a/dpnp/CMakeLists.txt +++ b/dpnp/CMakeLists.txt @@ -47,6 +47,7 @@ endfunction() build_dpnp_cython_ext_with_backend(dparray ${CMAKE_CURRENT_SOURCE_DIR}/dparray.pyx dpnp) add_subdirectory(backend) +add_subdirectory(backend/extensions/lapack) add_subdirectory(dpnp_algo) add_subdirectory(dpnp_utils) diff --git a/dpnp/backend/CMakeLists.txt b/dpnp/backend/CMakeLists.txt index 894d3b6e72d4..6fd13acd1128 100644 --- a/dpnp/backend/CMakeLists.txt +++ b/dpnp/backend/CMakeLists.txt @@ -109,7 +109,7 @@ add_library(dpnp_backend_library INTERFACE IMPORTED GLOBAL) target_include_directories(dpnp_backend_library BEFORE INTERFACE ${CMAKE_CURRENT_SOURCE_DIR}/include ${CMAKE_CURRENT_SOURCE_DIR}/src) target_link_libraries(dpnp_backend_library INTERFACE ${_trgt}) -if(DPNP_BACKEND_TESTS) +if (DPNP_BACKEND_TESTS) add_subdirectory(tests) endif() diff --git a/dpnp/backend/extensions/lapack/CMakeLists.txt b/dpnp/backend/extensions/lapack/CMakeLists.txt new file mode 100644 index 000000000000..a32adaa431ff --- /dev/null +++ b/dpnp/backend/extensions/lapack/CMakeLists.txt @@ -0,0 +1,75 @@ +# ***************************************************************************** +# Copyright (c) 2016-2023, Intel Corporation +# All rights reserved. +# +# Redistribution and use in source and binary forms, with or without +# modification, are permitted provided that the following conditions are met: +# - Redistributions of source code must retain the above copyright notice, +# this list of conditions and the following disclaimer. +# - 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. +# +# 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. +# ***************************************************************************** + + +set(python_module_name _lapack_impl) +pybind11_add_module(${python_module_name} MODULE + lapack_py.cpp + heevd.cpp + syevd.cpp +) + +if (WIN32) + if (${CMAKE_VERSION} VERSION_LESS "3.23") + # this is a work-around for target_link_options inserting option after -link option, cause + # linker to ignore it. + set(CMAKE_CXX_LINK_FLAGS "${CMAKE_CXX_LINK_FLAGS} -fsycl-device-code-split=per_kernel") + endif() +endif() + +set_target_properties(${python_module_name} PROPERTIES CMAKE_POSITION_INDEPENDENT_CODE ON) + +target_include_directories(${python_module_name} PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/../../include) +target_include_directories(${python_module_name} PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/../../src) + +target_include_directories(${python_module_name} PUBLIC ${Dpctl_INCLUDE_DIRS}) + +if (WIN32) + target_compile_options(${python_module_name} PRIVATE + /clang:-fno-approx-func + /clang:-fno-finite-math-only + ) +else() + target_compile_options(${python_module_name} PRIVATE + -fno-approx-func + -fno-finite-math-only + ) +endif() + +target_link_options(${python_module_name} PUBLIC -fsycl-device-code-split=per_kernel) +if (UNIX) + # this option is support on Linux only + target_link_options(${python_module_name} PUBLIC -fsycl-link-huge-device-code) +endif() + +if (DPNP_GENERATE_COVERAGE) + target_link_options(${python_module_name} PRIVATE -fprofile-instr-generate -fcoverage-mapping) +endif() + +target_link_libraries(${python_module_name} PUBLIC MKL::MKL_DPCPP) + +install(TARGETS ${python_module_name} + DESTINATION "dpnp/backend/extensions/lapack" +) diff --git a/dpnp/backend/extensions/lapack/heevd.cpp b/dpnp/backend/extensions/lapack/heevd.cpp new file mode 100644 index 000000000000..f873ee14d754 --- /dev/null +++ b/dpnp/backend/extensions/lapack/heevd.cpp @@ -0,0 +1,210 @@ +//***************************************************************************** +// Copyright (c) 2023, Intel Corporation +// All rights reserved. +// +// Redistribution and use in source and binary forms, with or without +// modification, are permitted provided that the following conditions are met: +// - Redistributions of source code must retain the above copyright notice, +// this list of conditions and the following disclaimer. +// - 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. +// +// 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. +//***************************************************************************** + + +#include + +#include "heevd.hpp" + +#include "dpnp_utils.hpp" + + +namespace dpnp +{ +namespace backend +{ +namespace ext +{ +namespace lapack +{ + +namespace mkl_lapack = oneapi::mkl::lapack; +namespace py = pybind11; + +template +static inline sycl::event call_heevd(sycl::queue exec_q, + const oneapi::mkl::job jobz, + const oneapi::mkl::uplo upper_lower, + const std::int64_t n, + T* a, + RealT* w, + std::vector &host_task_events, + const std::vector& depends) +{ + validate_type_for_device(exec_q); + validate_type_for_device(exec_q); + + const std::int64_t lda = std::max(1UL, n); + const std::int64_t scratchpad_size = mkl_lapack::heevd_scratchpad_size(exec_q, jobz, upper_lower, n, lda); + T* scratchpad = nullptr; + + std::stringstream error_msg; + std::int64_t info = 0; + + sycl::event heevd_event; + try + { + scratchpad = sycl::malloc_device(scratchpad_size, exec_q); + + heevd_event = mkl_lapack::heevd( + exec_q, + jobz, // 'jobz == job::vec' means eigenvalues and eigenvectors are computed. + upper_lower, // 'upper_lower == job::upper' means the upper triangular part of A, or the lower triangular otherwise + n, // The order of the matrix A (0 <= n) + a, // Pointer to A, size (lda, *), where the 2nd dimension, must be at least max(1, n) + // If 'jobz == job::vec', then on exit it will contain the eigenvectors of A + lda, // The leading dimension of a, must be at least max(1, n) + w, // Pointer to array of size at least n, it will contain the eigenvalues of A in ascending order + scratchpad, // Pointer to scratchpad memory to be used by MKL routine for storing intermediate results + scratchpad_size, + depends); + } + catch (mkl_lapack::exception const& e) + { + error_msg << "Unexpected MKL exception caught during heevd() call:\nreason: " << e.what() + << "\ninfo: " << e.info(); + info = e.info(); + } + catch (sycl::exception const& e) + { + error_msg << "Unexpected SYCL exception caught during heevd() call:\n" << e.what(); + info = -1; + } + + if (info != 0) // an unexected error occurs + { + if (scratchpad != nullptr) + { + sycl::free(scratchpad, exec_q); + } + throw std::runtime_error(error_msg.str()); + } + + sycl::event clean_up_event = exec_q.submit([&](sycl::handler& cgh) { + cgh.depends_on(heevd_event); + auto ctx = exec_q.get_context(); + cgh.host_task([ctx, scratchpad]() { sycl::free(scratchpad, ctx); }); + }); + host_task_events.push_back(clean_up_event); + return heevd_event; +} + +std::pair heevd(sycl::queue exec_q, + const std::int8_t jobz, + const std::int8_t upper_lower, + dpctl::tensor::usm_ndarray eig_vecs, + dpctl::tensor::usm_ndarray eig_vals, + const std::vector& depends) +{ + const int eig_vecs_nd = eig_vecs.get_ndim(); + const int eig_vals_nd = eig_vals.get_ndim(); + + if (eig_vecs_nd != 2) + { + throw py::value_error("Unexpected ndim=" + std::to_string(eig_vecs_nd) + + " of an output array with eigenvectors"); + } + else if (eig_vals_nd != 1) + { + throw py::value_error("Unexpected ndim=" + std::to_string(eig_vals_nd) + + " of an output array with eigenvalues"); + } + + const py::ssize_t* eig_vecs_shape = eig_vecs.get_shape_raw(); + const py::ssize_t* eig_vals_shape = eig_vals.get_shape_raw(); + + if (eig_vecs_shape[0] != eig_vecs_shape[1]) + { + throw py::value_error("Output array with eigenvectors with be square"); + } + else if (eig_vecs_shape[0] != eig_vals_shape[0]) + { + throw py::value_error("Eigenvectors and eigenvalues have different shapes"); + } + + size_t src_nelems(1); + + for (int i = 0; i < eig_vecs_nd; ++i) + { + src_nelems *= static_cast(eig_vecs_shape[i]); + } + + if (src_nelems == 0) + { + // nothing to do + return std::make_pair(sycl::event(), sycl::event()); + } + + // check compatibility of execution queue and allocation queue + if (!dpctl::utils::queues_are_compatible(exec_q, {eig_vecs, eig_vals})) + { + throw py::value_error("Execution queue is not compatible with allocation queues"); + } + + // check that arrays do not overlap, and concurrent access is safe. + // TODO: need to be exposed by DPCTL headers + // auto const &overlap = dpctl::tensor::overlap::MemoryOverlap(); + // if (overlap(eig_vecs, eig_vals)) + // { + // throw py::value_error("Arrays index overlapping segments of memory"); + // } + + int eig_vecs_typenum = eig_vecs.get_typenum(); + int eig_vals_typenum = eig_vals.get_typenum(); + auto const& dpctl_capi = dpctl::detail::dpctl_capi::get(); + + sycl::event heevd_ev; + std::vector host_task_events; + + const std::int64_t n = eig_vecs_shape[0]; + const oneapi::mkl::job jobz_val = static_cast(jobz); + const oneapi::mkl::uplo uplo_val = static_cast(upper_lower); + + if ((eig_vecs_typenum == dpctl_capi.UAR_CDOUBLE_) && (eig_vals_typenum == dpctl_capi.UAR_DOUBLE_)) + { + std::complex* a = reinterpret_cast*>(eig_vecs.get_data()); + double* w = reinterpret_cast(eig_vals.get_data()); + + heevd_ev = call_heevd(exec_q, jobz_val, uplo_val, n, a, w, host_task_events, depends); + } + else if ((eig_vecs_typenum == dpctl_capi.UAR_CFLOAT_) && (eig_vals_typenum == dpctl_capi.UAR_FLOAT_)) + { + std::complex* a = reinterpret_cast*>(eig_vecs.get_data()); + float* w = reinterpret_cast(eig_vals.get_data()); + + heevd_ev = call_heevd(exec_q, jobz_val, uplo_val, n, a, w, host_task_events, depends); + } + else + { + throw py::value_error("Unexpected types of either eigenvectors or eigenvalues"); + } + + sycl::event args_ev = dpctl::utils::keep_args_alive(exec_q, {eig_vecs, eig_vals}, host_task_events); + return std::make_pair(args_ev, heevd_ev); +} +} +} +} +} diff --git a/dpnp/backend/extensions/lapack/heevd.hpp b/dpnp/backend/extensions/lapack/heevd.hpp new file mode 100644 index 000000000000..93ce6fe560e1 --- /dev/null +++ b/dpnp/backend/extensions/lapack/heevd.hpp @@ -0,0 +1,51 @@ +//***************************************************************************** +// Copyright (c) 2023, Intel Corporation +// All rights reserved. +// +// Redistribution and use in source and binary forms, with or without +// modification, are permitted provided that the following conditions are met: +// - Redistributions of source code must retain the above copyright notice, +// this list of conditions and the following disclaimer. +// - 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. +// +// 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. +//***************************************************************************** + +#pragma once + +#include +#include + +#include + + +namespace dpnp +{ +namespace backend +{ +namespace ext +{ +namespace lapack +{ + extern std::pair heevd(sycl::queue exec_q, + const std::int8_t jobz, + const std::int8_t upper_lower, + dpctl::tensor::usm_ndarray eig_vecs, + dpctl::tensor::usm_ndarray eig_vals, + const std::vector& depends); +} +} +} +} diff --git a/dpnp/backend/extensions/lapack/lapack_py.cpp b/dpnp/backend/extensions/lapack/lapack_py.cpp new file mode 100644 index 000000000000..ea7506308032 --- /dev/null +++ b/dpnp/backend/extensions/lapack/lapack_py.cpp @@ -0,0 +1,57 @@ +//***************************************************************************** +// Copyright (c) 2023, Intel Corporation +// All rights reserved. +// +// Redistribution and use in source and binary forms, with or without +// modification, are permitted provided that the following conditions are met: +// - Redistributions of source code must retain the above copyright notice, +// this list of conditions and the following disclaimer. +// - 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. +// +// 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. +//***************************************************************************** +// +// This file defines functions of dpnp.backend._lapack_impl extensions +// +//***************************************************************************** + +#include +#include + +#include "heevd.hpp" +#include "syevd.hpp" + +namespace py = pybind11; + +PYBIND11_MODULE(_lapack_impl, m) +{ + m.def("_heevd", + &dpnp::backend::ext::lapack::heevd, + "Call `heevd` from OneMKL LAPACK library to return " + "the eigenvalues and eigenvectors of a complex Hermitian matrix", + py::arg("sycl_queue"), + py::arg("jobz"), py::arg("upper_lower"), + py::arg("eig_vecs"), py::arg("eig_vals"), + py::arg("depends") = py::list()); + + m.def("_syevd", + &dpnp::backend::ext::lapack::syevd, + "Call `syevd` from OneMKL LAPACK library to return " + "the eigenvalues and eigenvectors of a real symmetric matrix", + py::arg("sycl_queue"), + py::arg("jobz"), py::arg("upper_lower"), + py::arg("eig_vecs"), py::arg("eig_vals"), + py::arg("depends") = py::list()); +} diff --git a/dpnp/backend/extensions/lapack/syevd.cpp b/dpnp/backend/extensions/lapack/syevd.cpp new file mode 100644 index 000000000000..93be82d201d8 --- /dev/null +++ b/dpnp/backend/extensions/lapack/syevd.cpp @@ -0,0 +1,214 @@ +//***************************************************************************** +// Copyright (c) 2023, Intel Corporation +// All rights reserved. +// +// Redistribution and use in source and binary forms, with or without +// modification, are permitted provided that the following conditions are met: +// - Redistributions of source code must retain the above copyright notice, +// this list of conditions and the following disclaimer. +// - 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. +// +// 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. +//***************************************************************************** + + +#include + +#include "syevd.hpp" + +#include "dpnp_utils.hpp" + + +namespace dpnp +{ +namespace backend +{ +namespace ext +{ +namespace lapack +{ + +namespace mkl_lapack = oneapi::mkl::lapack; +namespace py = pybind11; + +template +static inline sycl::event call_syevd(sycl::queue exec_q, + const oneapi::mkl::job jobz, + const oneapi::mkl::uplo upper_lower, + const std::int64_t n, + T* a, + T* w, + std::vector &host_task_events, + const std::vector& depends) +{ + validate_type_for_device(exec_q); + + const std::int64_t lda = std::max(1UL, n); + const std::int64_t scratchpad_size = mkl_lapack::syevd_scratchpad_size(exec_q, jobz, upper_lower, n, lda); + T* scratchpad = nullptr; + + std::stringstream error_msg; + std::int64_t info = 0; + + sycl::event syevd_event; + try + { + scratchpad = sycl::malloc_device(scratchpad_size, exec_q); + + syevd_event = mkl_lapack::syevd( + exec_q, + jobz, // 'jobz == job::vec' means eigenvalues and eigenvectors are computed. + upper_lower, // 'upper_lower == job::upper' means the upper triangular part of A, or the lower triangular otherwise + n, // The order of the matrix A (0 <= n) + a, // Pointer to A, size (lda, *), where the 2nd dimension, must be at least max(1, n) + // If 'jobz == job::vec', then on exit it will contain the eigenvectors of A + lda, // The leading dimension of a, must be at least max(1, n) + w, // Pointer to array of size at least n, it will contain the eigenvalues of A in ascending order + scratchpad, // Pointer to scratchpad memory to be used by MKL routine for storing intermediate results + scratchpad_size, + depends); + } + catch (mkl_lapack::exception const& e) + { + error_msg << "Unexpected MKL exception caught during syevd() call:\nreason: " << e.what() + << "\ninfo: " << e.info(); + info = e.info(); + } + catch (sycl::exception const& e) + { + error_msg << "Unexpected SYCL exception caught during syevd() call:\n" << e.what(); + info = -1; + } + + if (info != 0) // an unexected error occurs + { + if (scratchpad != nullptr) + { + sycl::free(scratchpad, exec_q); + } + throw std::runtime_error(error_msg.str()); + } + + sycl::event clean_up_event = exec_q.submit([&](sycl::handler& cgh) { + cgh.depends_on(syevd_event); + auto ctx = exec_q.get_context(); + cgh.host_task([ctx, scratchpad]() { sycl::free(scratchpad, ctx); }); + }); + host_task_events.push_back(clean_up_event); + return syevd_event; +} + +std::pair syevd(sycl::queue exec_q, + const std::int8_t jobz, + const std::int8_t upper_lower, + dpctl::tensor::usm_ndarray eig_vecs, + dpctl::tensor::usm_ndarray eig_vals, + const std::vector& depends) +{ + const int eig_vecs_nd = eig_vecs.get_ndim(); + const int eig_vals_nd = eig_vals.get_ndim(); + + if (eig_vecs_nd != 2) + { + throw py::value_error("Unexpected ndim=" + std::to_string(eig_vecs_nd) + + " of an output array with eigenvectors"); + } + else if (eig_vals_nd != 1) + { + throw py::value_error("Unexpected ndim=" + std::to_string(eig_vals_nd) + + " of an output array with eigenvalues"); + } + + const py::ssize_t* eig_vecs_shape = eig_vecs.get_shape_raw(); + const py::ssize_t* eig_vals_shape = eig_vals.get_shape_raw(); + + if (eig_vecs_shape[0] != eig_vecs_shape[1]) + { + throw py::value_error("Output array with eigenvectors with be square"); + } + else if (eig_vecs_shape[0] != eig_vals_shape[0]) + { + throw py::value_error("Eigenvectors and eigenvalues have different shapes"); + } + + size_t src_nelems(1); + + for (int i = 0; i < eig_vecs_nd; ++i) + { + src_nelems *= static_cast(eig_vecs_shape[i]); + } + + if (src_nelems == 0) + { + // nothing to do + return std::make_pair(sycl::event(), sycl::event()); + } + + // check compatibility of execution queue and allocation queue + if (!dpctl::utils::queues_are_compatible(exec_q, {eig_vecs, eig_vals})) + { + throw py::value_error("Execution queue is not compatible with allocation queues"); + } + + // check that arrays do not overlap, and concurrent access is safe. + // TODO: need to be exposed by DPCTL headers + // auto const& overlap = dpctl::tensor::overlap::MemoryOverlap(); + // if (overlap(eig_vecs, eig_vals)) + // { + // throw py::value_error("Arrays index overlapping segments of memory"); + // } + + int eig_vecs_typenum = eig_vecs.get_typenum(); + int eig_vals_typenum = eig_vals.get_typenum(); + auto const& dpctl_capi = dpctl::detail::dpctl_capi::get(); + + sycl::event syevd_ev; + std::vector host_task_events; + + const std::int64_t n = eig_vecs_shape[0]; + const oneapi::mkl::job jobz_val = static_cast(jobz); + const oneapi::mkl::uplo uplo_val = static_cast(upper_lower); + + if (eig_vecs_typenum != eig_vals_typenum) + { + throw py::value_error("Types of eigenvectors and eigenvalues aare missmatched"); + } + else if (eig_vecs_typenum == dpctl_capi.UAR_DOUBLE_) + { + double* a = reinterpret_cast(eig_vecs.get_data()); + double* w = reinterpret_cast(eig_vals.get_data()); + + syevd_ev = call_syevd(exec_q, jobz_val, uplo_val, n, a, w, host_task_events, depends); + } + else if (eig_vecs_typenum == dpctl_capi.UAR_FLOAT_) + { + float* a = reinterpret_cast(eig_vecs.get_data()); + float* w = reinterpret_cast(eig_vals.get_data()); + + syevd_ev = call_syevd(exec_q, jobz_val, uplo_val, n, a, w, host_task_events, depends); + } + else + { + throw py::value_error("Unexpected types with num=" + std::to_string(eig_vecs_typenum) + + " for eigenvectors and eigenvalues"); + } + + sycl::event args_ev = dpctl::utils::keep_args_alive(exec_q, {eig_vecs, eig_vals}, host_task_events); + return std::make_pair(args_ev, syevd_ev); +} +} +} +} +} diff --git a/dpnp/backend/extensions/lapack/syevd.hpp b/dpnp/backend/extensions/lapack/syevd.hpp new file mode 100644 index 000000000000..14d167ec02a7 --- /dev/null +++ b/dpnp/backend/extensions/lapack/syevd.hpp @@ -0,0 +1,51 @@ +//***************************************************************************** +// Copyright (c) 2023, Intel Corporation +// All rights reserved. +// +// Redistribution and use in source and binary forms, with or without +// modification, are permitted provided that the following conditions are met: +// - Redistributions of source code must retain the above copyright notice, +// this list of conditions and the following disclaimer. +// - 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. +// +// 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. +//***************************************************************************** + +#pragma once + +#include +#include + +#include + + +namespace dpnp +{ +namespace backend +{ +namespace ext +{ +namespace lapack +{ + extern std::pair syevd(sycl::queue exec_q, + const std::int8_t jobz, + const std::int8_t upper_lower, + dpctl::tensor::usm_ndarray eig_vecs, + dpctl::tensor::usm_ndarray eig_vals, + const std::vector& depends = {}); +} +} +} +} diff --git a/dpnp/dpnp_iface.py b/dpnp/dpnp_iface.py index 6a5bcf239df2..f6aabbb1399d 100644 --- a/dpnp/dpnp_iface.py +++ b/dpnp/dpnp_iface.py @@ -67,7 +67,9 @@ "from_dlpack", "get_dpnp_descriptor", "get_include", - "get_normalized_queue_device" + "get_normalized_queue_device", + "get_usm_ndarray", + "is_supported_array_type" ] from dpnp import ( @@ -371,3 +373,44 @@ def get_normalized_queue_device(obj=None, if hasattr(dpt._device, 'normalize_queue_device'): return dpt._device.normalize_queue_device(sycl_queue=sycl_queue, device=device) return sycl_queue + + +def get_usm_ndarray(a): + """ + Return :class:`dpctl.tensor.usm_ndarray` from input array `a`. + + Parameters + ---------- + a : {dpnp_array, usm_ndarray} + Input array of supported type :class:`dpnp.ndarray` + or :class:`dpctl.tensor.usm_ndarray`. + + Returns + ------- + out : usm_ndarray + A dpctl USM ndarray of input array `a`. + + """ + + return a.get_array() if isinstance(a, dpnp_array) else a + + +def is_supported_array_type(a): + """ + Return ``True`` if an array of either type :class:`dpnp.ndarray` + or :class:`dpctl.tensor.usm_ndarray` type, ``False`` otherwise. + + Parameters + ---------- + a : array + An input array to check the type. + + Returns + ------- + out : bool + ``True`` if type of array `a` is supported array type, + ``False`` otherwise. + + """ + + return isinstance(a, (dpnp_array, dpt.usm_ndarray)) diff --git a/dpnp/linalg/dpnp_iface_linalg.py b/dpnp/linalg/dpnp_iface_linalg.py index 6e6f55db8f92..e2e962585786 100644 --- a/dpnp/linalg/dpnp_iface_linalg.py +++ b/dpnp/linalg/dpnp_iface_linalg.py @@ -41,6 +41,10 @@ import dpnp +from .dpnp_utils_linalg import ( + dpnp_eigh +) + import numpy from dpnp.dpnp_utils import * @@ -53,6 +57,7 @@ "cond", "det", "eig", + "eigh", "eigvals", "inv", "matrix_power", @@ -172,6 +177,68 @@ def eig(x1): return call_origin(numpy.linalg.eig, x1) +def eigh(a, UPLO='L'): + """ + Return the eigenvalues and eigenvectors of a complex Hermitian + (conjugate symmetric) or a real symmetric matrix. + + Returns two objects, a 1-D array containing the eigenvalues of `a`, and + a 2-D square array or matrix (depending on the input type) of the + corresponding eigenvectors (in columns). + + For full documentation refer to :obj:`numpy.linalg.eigh`. + + Returns + ------- + w : (..., M) dpnp.ndarray + The eigenvalues in ascending order, each repeated according to + its multiplicity. + v : (..., M, M) dpnp.ndarray + The column ``v[:, i]`` is the normalized eigenvector corresponding + to the eigenvalue ``w[i]``. + + Limitations + ----------- + Parameter `a` is supported as :class:`dpnp.ndarray` or :class:`dpctl.tensor.usm_ndarray`. + Input array data types are limited by supported DPNP :ref:`Data types`. + + See Also + -------- + :obj:`dpnp.eig` : eigenvalues and right eigenvectors for non-symmetric arrays. + :obj:`dpnp.eigvals` : eigenvalues of non-symmetric arrays. + + Examples + -------- + >>> import dpnp as dp + >>> a = dp.array([[1, -2j], [2j, 5]]) + >>> a + array([[ 1.+0.j, -0.-2.j], + [ 0.+2.j, 5.+0.j]]) + >>> w, v = dp.linalg.eigh(a) + >>> w; v + array([0.17157288, 5.82842712]), + array([[-0.92387953-0.j , -0.38268343+0.j ], # may vary + [ 0. +0.38268343j, 0. -0.92387953j]])) + + """ + + if UPLO not in ('L', 'U'): + raise ValueError("UPLO argument must be 'L' or 'U'") + + if not dpnp.is_supported_array_type(a): + raise TypeError("An array must be any of supported type, but got {}".format(type(a))) + + if a.ndim < 2: + raise ValueError("%d-dimensional array given. Array must be " + "at least two-dimensional" % a.ndim) + + m, n = a.shape[-2:] + if m != n: + raise ValueError("Last 2 dimensions of the array must be square") + + return dpnp_eigh(a, UPLO=UPLO) + + def eigvals(input): """ Compute the eigenvalues of a general matrix. diff --git a/dpnp/linalg/dpnp_utils_linalg.py b/dpnp/linalg/dpnp_utils_linalg.py new file mode 100644 index 000000000000..32a9ac7d5607 --- /dev/null +++ b/dpnp/linalg/dpnp_utils_linalg.py @@ -0,0 +1,134 @@ +# cython: language_level=3 +# distutils: language = c++ +# -*- coding: utf-8 -*- +# ***************************************************************************** +# Copyright (c) 2023, Intel Corporation +# All rights reserved. +# +# Redistribution and use in source and binary forms, with or without +# modification, are permitted provided that the following conditions are met: +# - Redistributions of source code must retain the above copyright notice, +# this list of conditions and the following disclaimer. +# - 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. +# +# 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. +# ***************************************************************************** + + +import dpnp +import dpnp.backend.extensions.lapack._lapack_impl as li + +import dpctl.tensor._tensor_impl as ti + +__all__ = [ + "dpnp_eigh" +] + +_jobz = {'N': 0, 'V': 1} +_upper_lower = {'U': 0, 'L': 1} + + +def dpnp_eigh(a, UPLO): + """ + Return the eigenvalues and eigenvectors of a complex Hermitian + (conjugate symmetric) or a real symmetric matrix. + + The main calculation is done by calling an extention function + for LAPACK library of OneMKL. Depending on input type of `a` array, + it will be either ``heevd`` (for complex types) or ``syevd`` (for others). + + """ + + a_usm_type = a.usm_type + a_sycl_queue = a.sycl_queue + a_order = 'C' if a.flags.c_contiguous else 'F' + a_usm_arr = dpnp.get_usm_ndarray(a) + + # 'V' means both eigenvectors and eigenvalues will be calculated + jobz = _jobz['V'] + uplo = _upper_lower[UPLO] + + # get resulting type of arrays with eigenvalues and eigenvectors + a_dtype = a.dtype + lapack_func = "_syevd" + if dpnp.issubdtype(a_dtype, dpnp.complexfloating): + lapack_func = "_heevd" + v_type = a_dtype + w_type = dpnp.float64 if a_dtype == dpnp.complex128 else dpnp.float32 + elif dpnp.issubdtype(a_dtype, dpnp.floating): + v_type = w_type = a_dtype + elif a_sycl_queue.sycl_device.has_aspect_fp64: + v_type = w_type = dpnp.float64 + else: + v_type = w_type = dpnp.float32 + + if a.ndim > 2: + w = dpnp.empty(a.shape[:-1], dtype=w_type, usm_type=a_usm_type, sycl_queue=a_sycl_queue) + + # need to loop over the 1st dimension to get eigenvalues and eigenvectors of 3d matrix A + op_count = a.shape[0] + if op_count == 0: + return w, dpnp.empty_like(a, dtype=v_type) + + eig_vecs = [None] * op_count + ht_copy_ev = [None] * op_count + ht_lapack_ev = [None] * op_count + for i in range(op_count): + # oneMKL LAPACK assumes fortran-like array as input, so + # allocate a memory with 'F' order for dpnp array of eigenvectors + eig_vecs[i] = dpnp.empty_like(a[i], order='F', dtype=v_type) + + # use DPCTL tensor function to fill the array of eigenvectors with content of input array + ht_copy_ev[i], copy_ev = ti._copy_usm_ndarray_into_usm_ndarray(src=a_usm_arr[i], dst=eig_vecs[i].get_array(), sycl_queue=a_sycl_queue) + + # call LAPACK extension function to get eigenvalues and eigenvectors of a portion of matrix A + ht_lapack_ev[i], _ = getattr(li, lapack_func)(a_sycl_queue, jobz, uplo, eig_vecs[i].get_array(), w[i].get_array(), depends=[copy_ev]) + + # TODO: remove once dpctl fix is available + ht_lapack_ev[i].wait() + + for i in range(op_count): + # ht_lapack_ev[i].wait() + ht_copy_ev[i].wait() + + # combine the list of eigenvectors into a single array + v = dpnp.array(eig_vecs, order=a_order) + return w, v + else: + # oneMKL LAPACK assumes fortran-like array as input, so + # allocate a memory with 'F' order for dpnp array of eigenvectors + v = dpnp.empty_like(a, order='F', dtype=v_type) + + # use DPCTL tensor function to fill the array of eigenvectors with content of input array + ht_copy_ev, copy_ev = ti._copy_usm_ndarray_into_usm_ndarray(src=a_usm_arr, dst=v.get_array(), sycl_queue=a_sycl_queue) + + # allocate a memory for dpnp array of eigenvalues + w = dpnp.empty(a.shape[:-1], dtype=w_type, usm_type=a_usm_type, sycl_queue=a_sycl_queue) + + # call LAPACK extension function to get eigenvalues and eigenvectors of matrix A + ht_lapack_ev, lapack_ev = getattr(li, lapack_func)(a_sycl_queue, jobz, uplo, v.get_array(), w.get_array(), depends=[copy_ev]) + + if a_order != 'F': + # need to align order of eigenvectors with one of input matrix A + out_v = dpnp.empty_like(v, order=a_order) + ht_copy_out_ev, _ = ti._copy_usm_ndarray_into_usm_ndarray(src=v.get_array(), dst=out_v.get_array(), sycl_queue=a_sycl_queue, depends=[lapack_ev]) + ht_copy_out_ev.wait() + else: + out_v = v + + ht_lapack_ev.wait() + ht_copy_ev.wait() + + return w, out_v diff --git a/tests/skipped_tests_gpu.tbl b/tests/skipped_tests_gpu.tbl index fee79df860fb..36b17d5edbc3 100644 --- a/tests/skipped_tests_gpu.tbl +++ b/tests/skipped_tests_gpu.tbl @@ -405,6 +405,8 @@ tests/third_party/cupy/core_tests/test_ndarray_reduction.py::TestArrayReduction: tests/third_party/cupy/core_tests/test_ndarray_reduction.py::TestArrayReduction::test_ptp_multiple_axes tests/third_party/cupy/core_tests/test_ndarray_reduction.py::TestArrayReduction::test_ptp_multiple_axes_keepdims tests/third_party/cupy/core_tests/test_ndarray_reduction.py::TestArrayReduction::test_ptp_nan +tests/third_party/cupy/core_tests/test_ndarray_reduction.py::TestArrayReduction::test_ptp_nan_imag +tests/third_party/cupy/core_tests/test_ndarray_reduction.py::TestArrayReduction::test_ptp_nan_real tests/third_party/cupy/core_tests/test_ndarray_reduction.py::TestCubReduction_param_0_{order='C', shape=(10,)}::test_cub_max tests/third_party/cupy/core_tests/test_ndarray_reduction.py::TestCubReduction_param_0_{order='C', shape=(10,)}::test_cub_min @@ -645,11 +647,14 @@ tests/third_party/cupy/indexing_tests/test_indexing.py::TestIndexing::test_take_ tests/third_party/cupy/indexing_tests/test_indexing.py::TestSelect::test_select tests/third_party/cupy/indexing_tests/test_indexing.py::TestSelect::test_select_1D_choicelist tests/third_party/cupy/indexing_tests/test_indexing.py::TestSelect::test_select_choicelist_condlist_broadcast +tests/third_party/cupy/indexing_tests/test_indexing.py::TestSelect::test_select_complex tests/third_party/cupy/indexing_tests/test_indexing.py::TestSelect::test_select_default +tests/third_party/cupy/indexing_tests/test_indexing.py::TestSelect::test_select_default_complex tests/third_party/cupy/indexing_tests/test_indexing.py::TestSelect::test_select_default_scalar tests/third_party/cupy/indexing_tests/test_indexing.py::TestSelect::test_select_empty_lists tests/third_party/cupy/indexing_tests/test_indexing.py::TestSelect::test_select_length_error tests/third_party/cupy/indexing_tests/test_indexing.py::TestSelect::test_select_odd_shaped_broadcastable +tests/third_party/cupy/indexing_tests/test_indexing.py::TestSelect::test_select_odd_shaped_broadcastable_complex tests/third_party/cupy/indexing_tests/test_indexing.py::TestSelect::test_select_odd_shaped_non_broadcastable tests/third_party/cupy/indexing_tests/test_indexing.py::TestSelect::test_select_type_error_condlist tests/third_party/cupy/indexing_tests/test_indexing.py::TestIndexing::test_diagonal @@ -992,9 +997,14 @@ tests/third_party/cupy/math_tests/test_rounding.py::TestRoundExtreme_param_5_{de tests/third_party/cupy/math_tests/test_rounding.py::TestRoundExtreme_param_5_{decimals=99}::test_round_small tests/third_party/cupy/math_tests/test_rounding.py::TestRoundExtreme_param_6_{decimals=100}::test_round_large tests/third_party/cupy/math_tests/test_rounding.py::TestRoundExtreme_param_6_{decimals=100}::test_round_small +tests/third_party/cupy/math_tests/test_rounding.py::TestRounding::test_around +tests/third_party/cupy/math_tests/test_rounding.py::TestRounding::test_ceil tests/third_party/cupy/math_tests/test_rounding.py::TestRounding::test_fix +tests/third_party/cupy/math_tests/test_rounding.py::TestRounding::test_floor tests/third_party/cupy/math_tests/test_rounding.py::TestRounding::test_rint tests/third_party/cupy/math_tests/test_rounding.py::TestRounding::test_rint_negative +tests/third_party/cupy/math_tests/test_rounding.py::TestRounding::test_round_ +tests/third_party/cupy/math_tests/test_rounding.py::TestRounding::test_trunc tests/third_party/cupy/math_tests/test_sumprod.py::TestSumprod::test_sum_all_transposed tests/third_party/cupy/math_tests/test_sumprod.py::TestSumprod::test_sum_all_transposed2 tests/third_party/cupy/math_tests/test_sumprod.py::TestSumprod::test_sum_axes @@ -1430,6 +1440,8 @@ tests/third_party/cupy/statistics_tests/test_histogram.py::TestHistogram::test_b tests/third_party/cupy/statistics_tests/test_histogram.py::TestHistogram::test_histogram tests/third_party/cupy/statistics_tests/test_histogram.py::TestHistogram::test_histogram_array_bins tests/third_party/cupy/statistics_tests/test_histogram.py::TestHistogram::test_histogram_bins_not_ordered +tests/third_party/cupy/statistics_tests/test_histogram.py::TestHistogram::test_histogram_complex_weights +tests/third_party/cupy/statistics_tests/test_histogram.py::TestHistogram::test_histogram_complex_weights_uneven_bins tests/third_party/cupy/statistics_tests/test_histogram.py::TestHistogram::test_histogram_density tests/third_party/cupy/statistics_tests/test_histogram.py::TestHistogram::test_histogram_empty tests/third_party/cupy/statistics_tests/test_histogram.py::TestHistogram::test_histogram_float_weights diff --git a/tests/test_linalg.py b/tests/test_linalg.py index d9784a41558f..d90ac8bf9c6c 100644 --- a/tests/test_linalg.py +++ b/tests/test_linalg.py @@ -1,5 +1,8 @@ import pytest -from .helper import get_all_dtypes +from .helper import ( + get_all_dtypes, + get_complex_dtypes +) import dpnp as inp @@ -33,12 +36,12 @@ def vvsort(val, vec, size, xp): unravel_imax = numpy.unravel_index(imax, val.shape) # swap elements in val array - temp = xp.array(val[unravel_i], dtype=vec.dtype, **val_kwargs) + temp = xp.array(val[unravel_i], dtype=val.dtype, **val_kwargs) val[unravel_i] = val[unravel_imax] val[unravel_imax] = temp # swap corresponding columns in vec matrix - temp = xp.array(vec[:, i], dtype=val.dtype, **vec_kwargs) + temp = xp.array(vec[:, i], dtype=vec.dtype, **vec_kwargs) vec[:, i] = vec[:, imax] vec[:, imax] = temp @@ -126,11 +129,49 @@ def test_eig_arange(type, size): assert (dpnp_vec.dtype == np_vec.dtype) assert (dpnp_val.shape == np_val.shape) assert (dpnp_vec.shape == np_vec.shape) + assert (dpnp_val.usm_type == dpnp_symm.usm_type) + assert (dpnp_vec.usm_type == dpnp_symm.usm_type) assert_allclose(dpnp_val, np_val, rtol=1e-05, atol=1e-05) assert_allclose(dpnp_vec, np_vec, rtol=1e-05, atol=1e-05) +@pytest.mark.usefixtures("allow_fall_back_on_numpy") +@pytest.mark.parametrize("type", get_all_dtypes(no_bool=True, no_none=True)) +@pytest.mark.parametrize("size", [2, 4, 8]) +def test_eigh_arange(type, size): + a = numpy.arange(size * size, dtype=type).reshape((size, size)) + symm_orig = numpy.tril(a) + numpy.tril(a, -1).T + numpy.diag(numpy.full((size,), size * size, dtype=type)) + symm = symm_orig + dpnp_symm_orig = inp.array(symm) + dpnp_symm = dpnp_symm_orig + + dpnp_val, dpnp_vec = inp.linalg.eigh(dpnp_symm) + np_val, np_vec = numpy.linalg.eigh(symm) + + # DPNP sort val/vec by abs value + vvsort(dpnp_val, dpnp_vec, size, inp) + + # NP sort val/vec by abs value + vvsort(np_val, np_vec, size, numpy) + + # NP change sign of vectors + for i in range(np_vec.shape[1]): + if (np_vec[0, i] * dpnp_vec[0, i]).asnumpy() < 0: + np_vec[:, i] = -np_vec[:, i] + + assert_array_equal(symm_orig, symm) + assert_array_equal(dpnp_symm_orig, dpnp_symm) + + assert (dpnp_val.shape == np_val.shape) + assert (dpnp_vec.shape == np_vec.shape) + assert (dpnp_val.usm_type == dpnp_symm.usm_type) + assert (dpnp_vec.usm_type == dpnp_symm.usm_type) + + assert_allclose(dpnp_val, np_val, rtol=1e-05, atol=1e-04) + assert_allclose(dpnp_vec, np_vec, rtol=1e-05, atol=1e-04) + + @pytest.mark.parametrize("type", get_all_dtypes(no_bool=True, no_complex=True)) def test_eigvals(type): if dpctl.get_current_device_type() != dpctl.device_type.gpu: diff --git a/tests/test_sycl_queue.py b/tests/test_sycl_queue.py index fcea0d82eb86..e5e53646e1a1 100644 --- a/tests/test_sycl_queue.py +++ b/tests/test_sycl_queue.py @@ -769,6 +769,38 @@ def test_eig(device): assert_sycl_queue_equal(dpnp_vec_queue, expected_queue) +@pytest.mark.usefixtures("allow_fall_back_on_numpy") +@pytest.mark.parametrize("device", + valid_devices, + ids=[device.filter_string for device in valid_devices]) +def test_eigh(device): + size = 4 + a = numpy.arange(size * size, dtype=numpy.float64).reshape((size, size)) + symm_orig = numpy.tril(a) + numpy.tril(a, -1).T + numpy.diag(numpy.full((size,), size * size, dtype=numpy.float64)) + numpy_data = symm_orig + dpnp_symm_orig = dpnp.array(numpy_data, device=device) + dpnp_data = dpnp_symm_orig + + dpnp_val, dpnp_vec = dpnp.linalg.eigh(dpnp_data) + numpy_val, numpy_vec = numpy.linalg.eigh(numpy_data) + + assert_allclose(dpnp_val, numpy_val, rtol=1e-05, atol=1e-05) + assert_allclose(dpnp_vec, numpy_vec, rtol=1e-05, atol=1e-05) + + assert (dpnp_val.dtype == numpy_val.dtype) + assert (dpnp_vec.dtype == numpy_vec.dtype) + assert (dpnp_val.shape == numpy_val.shape) + assert (dpnp_vec.shape == numpy_vec.shape) + + expected_queue = dpnp_data.get_array().sycl_queue + dpnp_val_queue = dpnp_val.get_array().sycl_queue + dpnp_vec_queue = dpnp_vec.get_array().sycl_queue + + # compare queue and device + assert_sycl_queue_equal(dpnp_val_queue, expected_queue) + assert_sycl_queue_equal(dpnp_vec_queue, expected_queue) + + @pytest.mark.parametrize("device", valid_devices, ids=[device.filter_string for device in valid_devices]) diff --git a/tests/third_party/cupy/creation_tests/test_ranges.py b/tests/third_party/cupy/creation_tests/test_ranges.py index ac94297354f0..11e1d7f96048 100644 --- a/tests/third_party/cupy/creation_tests/test_ranges.py +++ b/tests/third_party/cupy/creation_tests/test_ranges.py @@ -192,7 +192,7 @@ def test_linspace_array_start_stop_axis1(self, xp, dtype_range, dtype_out): @testing.with_requires('numpy>=1.16') @testing.for_complex_dtypes() - @testing.numpy_cupy_array_equal() + @testing.numpy_cupy_allclose() def test_linspace_complex_start_stop(self, xp, dtype): start = xp.array([0, 120], dtype=dtype) stop = xp.array([100, 0], dtype=dtype) diff --git a/tests/third_party/cupy/linalg_tests/test_eigenvalue.py b/tests/third_party/cupy/linalg_tests/test_eigenvalue.py new file mode 100644 index 000000000000..fe577e32b285 --- /dev/null +++ b/tests/third_party/cupy/linalg_tests/test_eigenvalue.py @@ -0,0 +1,199 @@ +import unittest + +import numpy +import pytest + +import dpnp as cupy +from tests.third_party.cupy import testing + + +def _get_hermitian(xp, a, UPLO): + # TODO: fix this, currently dpnp.transpose() doesn't support complex types + # and no dpnp_array.swapaxes() + a = _wrap_as_numpy_array(xp, a) + _xp = numpy + + if UPLO == 'U': + _a = _xp.triu(a) + _xp.triu(a, k=1).swapaxes(-2, -1).conj() + else: + _a = _xp.tril(a) + _xp.tril(a, k=-1).swapaxes(-2, -1).conj() + return xp.array(_a) + +# TODO: remove once all required functionality is supported +def _wrap_as_numpy_array(xp, a): + return a.asnumpy() if xp is cupy else a + +@testing.parameterize(*testing.product({ + 'UPLO': ['U', 'L'], +})) +class TestEigenvalue(unittest.TestCase): + + @testing.for_all_dtypes() + @testing.numpy_cupy_allclose(rtol=1e-3, atol=1e-4, contiguous_check=False) + def test_eigh(self, xp, dtype): + if xp == numpy and dtype == numpy.float16: + # NumPy's eigh does not support float16 + _dtype = 'f' + else: + _dtype = dtype + if numpy.dtype(_dtype).kind == 'c': + a = xp.array([[1, 2j, 3], [4j, 5, 6j], [7, 8j, 9]], _dtype) + else: + a = xp.array([[1, 0, 3], [0, 5, 0], [7, 0, 9]], _dtype) + w, v = xp.linalg.eigh(a, UPLO=self.UPLO) + + # Changed the verification method to check if Av and vw match, since + # the eigenvectors of eigh() with CUDA 11.6 are mathematically correct + # but may not match NumPy. + A = _get_hermitian(xp, a, self.UPLO) + if _dtype == numpy.float16: + tol = 1e-3 + else: + tol = 1e-5 + + # TODO: remove _wrap_as_numpy_array() once @ support complex types + testing.assert_allclose(_wrap_as_numpy_array(xp, A) @ _wrap_as_numpy_array(xp, v), + _wrap_as_numpy_array(xp, v) @ numpy.diag(_wrap_as_numpy_array(xp, w)), + atol=tol, rtol=tol) + + # Check if v @ vt is an identity matrix + testing.assert_allclose(_wrap_as_numpy_array(xp, v) @ _wrap_as_numpy_array(xp, v).swapaxes(-2, -1).conj(), + numpy.identity(_wrap_as_numpy_array(xp, A).shape[-1], _dtype), atol=tol, + rtol=tol) + if xp == numpy and dtype == numpy.float16: + w = w.astype('e') + return w + + @testing.for_all_dtypes(no_bool=True, no_float16=True, no_complex=True) + @testing.numpy_cupy_allclose(rtol=1e-3, atol=1e-4, contiguous_check=False) + def test_eigh_batched(self, xp, dtype): + a = xp.array([[[1, 0, 3], [0, 5, 0], [7, 0, 9]], + [[3, 0, 3], [0, 7, 0], [7, 0, 11]]], dtype) + w, v = xp.linalg.eigh(a, UPLO=self.UPLO) + + # NumPy, cuSOLVER, rocSOLVER all sort in ascending order, + # so w's should be directly comparable. However, both cuSOLVER + # and rocSOLVER pick a different convention for constructing + # eigenvectors, so v's are not directly comparible and we verify + # them through the eigen equation A*v=w*v. + A = _get_hermitian(xp, a, self.UPLO) + for i in range(a.shape[0]): + testing.assert_allclose( + A[i].dot(v[i]), w[i]*v[i], rtol=1e-5, atol=1e-5) + return w + + @testing.for_complex_dtypes() + @testing.numpy_cupy_allclose(rtol=1e-3, atol=1e-4, contiguous_check=False) + def test_eigh_complex_batched(self, xp, dtype): + print() + a = xp.array([[[1, 2j, 3], [4j, 5, 6j], [7, 8j, 9]], + [[0, 2j, 3], [4j, 4, 6j], [7, 8j, 8]]], dtype) + w, v = xp.linalg.eigh(a, UPLO=self.UPLO) + + # NumPy, cuSOLVER, rocSOLVER all sort in ascending order, + # so w's should be directly comparable. However, both cuSOLVER + # and rocSOLVER pick a different convention for constructing + # eigenvectors, so v's are not directly comparible and we verify + # them through the eigen equation A*v=w*v. + A = _get_hermitian(xp, a, self.UPLO) + + # TODO: remove _wrap_as_numpy_array() once dpnp.dot() support complex types + A = _wrap_as_numpy_array(xp, A) + v = _wrap_as_numpy_array(xp, v) + w = _wrap_as_numpy_array(xp, w) + + for i in range(a.shape[0]): + testing.assert_allclose( + A[i].dot(v[i]), w[i]*v[i], rtol=1e-5, atol=1e-5) + return w + + @pytest.mark.skip("No support of dpnp.eigvalsh()") + @testing.for_all_dtypes(no_float16=True, no_complex=True) + @testing.numpy_cupy_allclose(rtol=1e-3, atol=1e-4) + def test_eigvalsh(self, xp, dtype): + a = xp.array([[1, 0, 3], [0, 5, 0], [7, 0, 9]], dtype) + w = xp.linalg.eigvalsh(a, UPLO=self.UPLO) + # NumPy, cuSOLVER, rocSOLVER all sort in ascending order, + # so they should be directly comparable + return w + + @pytest.mark.skip("No support of dpnp.eigvalsh()") + @testing.for_all_dtypes(no_float16=True, no_complex=True) + @testing.numpy_cupy_allclose(rtol=1e-3, atol=1e-4) + def test_eigvalsh_batched(self, xp, dtype): + a = xp.array([[[1, 0, 3], [0, 5, 0], [7, 0, 9]], + [[3, 0, 3], [0, 7, 0], [7, 0, 11]]], dtype) + w = xp.linalg.eigvalsh(a, UPLO=self.UPLO) + # NumPy, cuSOLVER, rocSOLVER all sort in ascending order, + # so they should be directly comparable + return w + + @pytest.mark.skip("No support of dpnp.eigvalsh()") + @testing.for_complex_dtypes() + @testing.numpy_cupy_allclose(rtol=1e-3, atol=1e-4) + def test_eigvalsh_complex(self, xp, dtype): + a = xp.array([[1, 2j, 3], [4j, 5, 6j], [7, 8j, 9]], dtype) + w = xp.linalg.eigvalsh(a, UPLO=self.UPLO) + # NumPy, cuSOLVER, rocSOLVER all sort in ascending order, + # so they should be directly comparable + return w + + @pytest.mark.skip("No support of dpnp.eigvalsh()") + @testing.for_complex_dtypes() + @testing.numpy_cupy_allclose(rtol=1e-3, atol=1e-4) + def test_eigvalsh_complex_batched(self, xp, dtype): + a = xp.array([[[1, 2j, 3], [4j, 5, 6j], [7, 8j, 9]], + [[0, 2j, 3], [4j, 4, 6j], [7, 8j, 8]]], dtype) + w = xp.linalg.eigvalsh(a, UPLO=self.UPLO) + # NumPy, cuSOLVER, rocSOLVER all sort in ascending order, + # so they should be directly comparable + return w + + +@testing.parameterize(*testing.product({ + 'UPLO': ['U', 'L'], + 'shape': [(0, 0), + (2, 0, 0), + (0, 3, 3)] +})) +class TestEigenvalueEmpty(unittest.TestCase): + + @testing.for_dtypes('ifdFD') + @testing.numpy_cupy_allclose() + def test_eigh(self, xp, dtype): + a = xp.empty(self.shape, dtype=dtype) + assert a.size == 0 + return xp.linalg.eigh(a, UPLO=self.UPLO) + + @pytest.mark.skip("No support of dpnp.eigvalsh()") + @testing.for_dtypes('ifdFD') + @testing.numpy_cupy_allclose() + def test_eigvalsh(self, xp, dtype): + a = xp.empty(self.shape, dtype) + assert a.size == 0 + return xp.linalg.eigvalsh(a, UPLO=self.UPLO) + + +@testing.parameterize(*testing.product({ + 'UPLO': ['U', 'L'], + 'shape': [(), + (3,), + (2, 3), + (4, 0), + (2, 2, 3), + (0, 2, 3)] +})) +class TestEigenvalueInvalid(unittest.TestCase): + + def test_eigh_shape_error(self): + for xp in (numpy, cupy): + a = xp.zeros(self.shape) + with pytest.raises((numpy.linalg.LinAlgError, ValueError)): + xp.linalg.eigh(a, self.UPLO) + + @pytest.mark.skip("No support of dpnp.eigvalsh()") + def test_eigvalsh_shape_error(self): + for xp in (numpy, cupy): + a = xp.zeros(self.shape) + with pytest.raises((numpy.linalg.LinAlgError, ValueError)): + xp.linalg.eigvalsh(a, self.UPLO) diff --git a/tests/third_party/cupy/testing/helper.py b/tests/third_party/cupy/testing/helper.py index 6331309820d2..af8f6e545b29 100644 --- a/tests/third_party/cupy/testing/helper.py +++ b/tests/third_party/cupy/testing/helper.py @@ -200,6 +200,16 @@ def _contains_signed_and_unsigned(kw): any(d in vs for d in _float_dtypes + _signed_dtypes) +def _wraps_partial(wrapped, *names): + # Only `wrapped` function have args of `names`. + def decorator(impl): + impl = functools.wraps(wrapped)(impl) + impl.__signature__ = inspect.signature( + functools.partial(wrapped, **{name: None for name in names})) + return impl + return decorator + + def _make_decorator(check_func, name, type_check, accept_error, sp_name=None, scipy_name=None): assert isinstance(name, str) @@ -640,16 +650,16 @@ def for_dtypes(dtypes, name='dtype'): argument. """ def decorator(impl): - @functools.wraps(impl) - def test_func(self, *args, **kw): + @_wraps_partial(impl, name) + def test_func(*args, **kw): for dtype in dtypes: try: kw[name] = numpy.dtype(dtype).type - impl(self, *args, **kw) + impl(*args, **kw) except unittest.SkipTest as e: - pass # print(f"Function decorator(): skipped: name={name} dtype={dtype} error={e}") + print('skipped: {} = {} ({})'.format(name, dtype, e)) except Exception: - # print(f"Function decorator(): name={name} dtype={dtype}") + print(name, 'is', dtype) raise return test_func @@ -661,8 +671,14 @@ def _get_supported_float_dtypes(): else: return (numpy.float32,) +def _get_supported_complex_dtypes(): + if select_default_device().has_aspect_fp64: + return (numpy.complex128, numpy.complex64) + else: + return (numpy.complex64,) + -_complex_dtypes = () +_complex_dtypes = _get_supported_complex_dtypes() _regular_float_dtypes = _get_supported_float_dtypes() _float_dtypes = _regular_float_dtypes _signed_dtypes = () From 000f9d63a4496970df89006da93382ef9521d69f Mon Sep 17 00:00:00 2001 From: Anton Volkov Date: Thu, 20 Apr 2023 15:15:58 -0500 Subject: [PATCH 2/3] Applying the review comments --- .gitignore | 1 - dpnp/backend/extensions/lapack/heevd.cpp | 27 +++++++++++++++++------- dpnp/backend/extensions/lapack/syevd.cpp | 27 +++++++++++++++++------- dpnp/linalg/dpnp_utils_linalg.py | 5 +---- 4 files changed, 39 insertions(+), 21 deletions(-) diff --git a/.gitignore b/.gitignore index ea9f2cba333d..8beb38f1efd6 100644 --- a/.gitignore +++ b/.gitignore @@ -1,6 +1,5 @@ # CMake build and local install directory _skbuild -build build_cython dpnp.egg-info diff --git a/dpnp/backend/extensions/lapack/heevd.cpp b/dpnp/backend/extensions/lapack/heevd.cpp index f873ee14d754..8c943646ff0a 100644 --- a/dpnp/backend/extensions/lapack/heevd.cpp +++ b/dpnp/backend/extensions/lapack/heevd.cpp @@ -44,14 +44,14 @@ namespace mkl_lapack = oneapi::mkl::lapack; namespace py = pybind11; template -static inline sycl::event call_heevd(sycl::queue exec_q, - const oneapi::mkl::job jobz, - const oneapi::mkl::uplo upper_lower, - const std::int64_t n, - T* a, - RealT* w, - std::vector &host_task_events, - const std::vector& depends) +static sycl::event call_heevd(sycl::queue exec_q, + const oneapi::mkl::job jobz, + const oneapi::mkl::uplo upper_lower, + const std::int64_t n, + T* a, + RealT* w, + std::vector& host_task_events, + const std::vector& depends) { validate_type_for_device(exec_q); validate_type_for_device(exec_q); @@ -171,6 +171,17 @@ std::pair heevd(sycl::queue exec_q, // throw py::value_error("Arrays index overlapping segments of memory"); // } + bool is_eig_vecs_f_contig = eig_vecs.is_f_contiguous(); + bool is_eig_vals_c_contig = eig_vals.is_c_contiguous(); + if (!is_eig_vecs_f_contig) + { + throw py::value_error("An array with input matrix / ouput eigenvectors must be F-contiguous"); + } + else if (!is_eig_vals_c_contig) + { + throw py::value_error("An array with output eigenvalues must be C-contiguous"); + } + int eig_vecs_typenum = eig_vecs.get_typenum(); int eig_vals_typenum = eig_vals.get_typenum(); auto const& dpctl_capi = dpctl::detail::dpctl_capi::get(); diff --git a/dpnp/backend/extensions/lapack/syevd.cpp b/dpnp/backend/extensions/lapack/syevd.cpp index 93be82d201d8..a4dded7543ab 100644 --- a/dpnp/backend/extensions/lapack/syevd.cpp +++ b/dpnp/backend/extensions/lapack/syevd.cpp @@ -44,14 +44,14 @@ namespace mkl_lapack = oneapi::mkl::lapack; namespace py = pybind11; template -static inline sycl::event call_syevd(sycl::queue exec_q, - const oneapi::mkl::job jobz, - const oneapi::mkl::uplo upper_lower, - const std::int64_t n, - T* a, - T* w, - std::vector &host_task_events, - const std::vector& depends) +static sycl::event call_syevd(sycl::queue exec_q, + const oneapi::mkl::job jobz, + const oneapi::mkl::uplo upper_lower, + const std::int64_t n, + T* a, + T* w, + std::vector& host_task_events, + const std::vector& depends) { validate_type_for_device(exec_q); @@ -170,6 +170,17 @@ std::pair syevd(sycl::queue exec_q, // throw py::value_error("Arrays index overlapping segments of memory"); // } + bool is_eig_vecs_f_contig = eig_vecs.is_f_contiguous(); + bool is_eig_vals_c_contig = eig_vals.is_c_contiguous(); + if (!is_eig_vecs_f_contig) + { + throw py::value_error("An array with input matrix / ouput eigenvectors must be F-contiguous"); + } + else if (!is_eig_vals_c_contig) + { + throw py::value_error("An array with output eigenvalues must be C-contiguous"); + } + int eig_vecs_typenum = eig_vecs.get_typenum(); int eig_vals_typenum = eig_vals.get_typenum(); auto const& dpctl_capi = dpctl::detail::dpctl_capi::get(); diff --git a/dpnp/linalg/dpnp_utils_linalg.py b/dpnp/linalg/dpnp_utils_linalg.py index 32a9ac7d5607..b7218b75d817 100644 --- a/dpnp/linalg/dpnp_utils_linalg.py +++ b/dpnp/linalg/dpnp_utils_linalg.py @@ -96,11 +96,8 @@ def dpnp_eigh(a, UPLO): # call LAPACK extension function to get eigenvalues and eigenvectors of a portion of matrix A ht_lapack_ev[i], _ = getattr(li, lapack_func)(a_sycl_queue, jobz, uplo, eig_vecs[i].get_array(), w[i].get_array(), depends=[copy_ev]) - # TODO: remove once dpctl fix is available - ht_lapack_ev[i].wait() - for i in range(op_count): - # ht_lapack_ev[i].wait() + ht_lapack_ev[i].wait() ht_copy_ev[i].wait() # combine the list of eigenvectors into a single array From 3a803f466587d0c2fa57e198d6a2bd44709ae9b3 Mon Sep 17 00:00:00 2001 From: Anton Volkov Date: Sat, 29 Apr 2023 14:25:53 +0200 Subject: [PATCH 3/3] Added array type check in dpnp.get_usm_ndarray() --- dpnp/dpnp_iface.py | 11 ++++++++++- 1 file changed, 10 insertions(+), 1 deletion(-) diff --git a/dpnp/dpnp_iface.py b/dpnp/dpnp_iface.py index f6aabbb1399d..ce3c540539d6 100644 --- a/dpnp/dpnp_iface.py +++ b/dpnp/dpnp_iface.py @@ -390,9 +390,18 @@ def get_usm_ndarray(a): out : usm_ndarray A dpctl USM ndarray of input array `a`. + Raises + ------ + TypeError + If input parameter `a` is of unsupported array type. + """ - return a.get_array() if isinstance(a, dpnp_array) else a + if isinstance(a, dpnp_array): + return a.get_array() + if isinstance(a, dpt.usm_ndarray): + return a + raise TypeError("An array must be any of supported type, but got {}".format(type(a))) def is_supported_array_type(a):