diff --git a/.github/workflows/conda-package.yml b/.github/workflows/conda-package.yml index b633b4bdc88..e50804468d8 100644 --- a/.github/workflows/conda-package.yml +++ b/.github/workflows/conda-package.yml @@ -30,6 +30,7 @@ env: test_umath.py test_usm_type.py third_party/cupy/core_tests + third_party/cupy/linalg_tests/test_decomposition.py third_party/cupy/linalg_tests/test_norms.py third_party/cupy/linalg_tests/test_product.py third_party/cupy/linalg_tests/test_solve.py diff --git a/dpnp/backend/extensions/lapack/CMakeLists.txt b/dpnp/backend/extensions/lapack/CMakeLists.txt index d8c902a1cf5..341845aa2db 100644 --- a/dpnp/backend/extensions/lapack/CMakeLists.txt +++ b/dpnp/backend/extensions/lapack/CMakeLists.txt @@ -31,6 +31,8 @@ set(_module_src ${CMAKE_CURRENT_SOURCE_DIR}/getrf.cpp ${CMAKE_CURRENT_SOURCE_DIR}/getrf_batch.cpp ${CMAKE_CURRENT_SOURCE_DIR}/heevd.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/potrf.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/potrf_batch.cpp ${CMAKE_CURRENT_SOURCE_DIR}/syevd.cpp ) diff --git a/dpnp/backend/extensions/lapack/getrf.cpp b/dpnp/backend/extensions/lapack/getrf.cpp index faebb9f7d42..f97d395bcd6 100644 --- a/dpnp/backend/extensions/lapack/getrf.cpp +++ b/dpnp/backend/extensions/lapack/getrf.cpp @@ -1,5 +1,5 @@ //***************************************************************************** -// Copyright (c) 2023, Intel Corporation +// Copyright (c) 2024, Intel Corporation // All rights reserved. // // Redistribution and use in source and binary forms, with or without diff --git a/dpnp/backend/extensions/lapack/getrf.hpp b/dpnp/backend/extensions/lapack/getrf.hpp index d8ea425f680..fee9b209426 100644 --- a/dpnp/backend/extensions/lapack/getrf.hpp +++ b/dpnp/backend/extensions/lapack/getrf.hpp @@ -1,5 +1,5 @@ //***************************************************************************** -// Copyright (c) 2023, Intel Corporation +// Copyright (c) 2024, Intel Corporation // All rights reserved. // // Redistribution and use in source and binary forms, with or without diff --git a/dpnp/backend/extensions/lapack/getrf_batch.cpp b/dpnp/backend/extensions/lapack/getrf_batch.cpp index 0944d7437ae..76977bf6628 100644 --- a/dpnp/backend/extensions/lapack/getrf_batch.cpp +++ b/dpnp/backend/extensions/lapack/getrf_batch.cpp @@ -1,5 +1,5 @@ //***************************************************************************** -// Copyright (c) 2023, Intel Corporation +// Copyright (c) 2024, Intel Corporation // All rights reserved. // // Redistribution and use in source and binary forms, with or without diff --git a/dpnp/backend/extensions/lapack/lapack_py.cpp b/dpnp/backend/extensions/lapack/lapack_py.cpp index e45cb685a1b..f97a0dc4a43 100644 --- a/dpnp/backend/extensions/lapack/lapack_py.cpp +++ b/dpnp/backend/extensions/lapack/lapack_py.cpp @@ -34,6 +34,7 @@ #include "getrf.hpp" #include "heevd.hpp" #include "linalg_exceptions.hpp" +#include "potrf.hpp" #include "syevd.hpp" namespace lapack_ext = dpnp::backend::ext::lapack; @@ -45,6 +46,8 @@ void init_dispatch_vectors(void) lapack_ext::init_gesv_dispatch_vector(); lapack_ext::init_getrf_batch_dispatch_vector(); lapack_ext::init_getrf_dispatch_vector(); + lapack_ext::init_potrf_batch_dispatch_vector(); + lapack_ext::init_potrf_dispatch_vector(); lapack_ext::init_syevd_dispatch_vector(); } @@ -92,6 +95,20 @@ PYBIND11_MODULE(_lapack_impl, m) py::arg("eig_vecs"), py::arg("eig_vals"), py::arg("depends") = py::list()); + m.def("_potrf", &lapack_ext::potrf, + "Call `potrf` from OneMKL LAPACK library to return " + "the Cholesky factorization of a symmetric positive-definite matrix", + py::arg("sycl_queue"), py::arg("a_array"), py::arg("upper_lower"), + py::arg("depends") = py::list()); + + m.def("_potrf_batch", &lapack_ext::potrf_batch, + "Call `potrf_batch` from OneMKL LAPACK library to return " + "the Cholesky factorization of a batch of symmetric " + "positive-definite matrix", + py::arg("sycl_queue"), py::arg("a_array"), py::arg("upper_lower"), + py::arg("n"), py::arg("stride_a"), py::arg("batch_size"), + py::arg("depends") = py::list()); + m.def("_syevd", &lapack_ext::syevd, "Call `syevd` from OneMKL LAPACK library to return " "the eigenvalues and eigenvectors of a real symmetric matrix", diff --git a/dpnp/backend/extensions/lapack/potrf.cpp b/dpnp/backend/extensions/lapack/potrf.cpp new file mode 100644 index 00000000000..610a629a9eb --- /dev/null +++ b/dpnp/backend/extensions/lapack/potrf.cpp @@ -0,0 +1,221 @@ +//***************************************************************************** +// Copyright (c) 2024, 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 + +// dpctl tensor headers +#include "utils/memory_overlap.hpp" +#include "utils/type_utils.hpp" + +#include "linalg_exceptions.hpp" +#include "potrf.hpp" +#include "types_matrix.hpp" + +#include "dpnp_utils.hpp" + +namespace dpnp +{ +namespace backend +{ +namespace ext +{ +namespace lapack +{ +namespace mkl_lapack = oneapi::mkl::lapack; +namespace py = pybind11; +namespace type_utils = dpctl::tensor::type_utils; + +typedef sycl::event (*potrf_impl_fn_ptr_t)(sycl::queue, + const oneapi::mkl::uplo, + const std::int64_t, + char *, + std::int64_t, + std::vector &, + const std::vector &); + +static potrf_impl_fn_ptr_t potrf_dispatch_vector[dpctl_td_ns::num_types]; + +template +static sycl::event potrf_impl(sycl::queue exec_q, + const oneapi::mkl::uplo upper_lower, + const std::int64_t n, + char *in_a, + std::int64_t lda, + std::vector &host_task_events, + const std::vector &depends) +{ + type_utils::validate_type_for_device(exec_q); + + T *a = reinterpret_cast(in_a); + + const std::int64_t scratchpad_size = + mkl_lapack::potrf_scratchpad_size(exec_q, upper_lower, n, lda); + T *scratchpad = nullptr; + + std::stringstream error_msg; + std::int64_t info = 0; + bool is_exception_caught = false; + + sycl::event potrf_event; + try { + scratchpad = sycl::malloc_device(scratchpad_size, exec_q); + + potrf_event = mkl_lapack::potrf( + exec_q, + upper_lower, // An enumeration value of type oneapi::mkl::uplo: + // oneapi::mkl::uplo::upper for the upper triangular + // part; oneapi::mkl::uplo::lower for the lower + // triangular part. + n, // Order of the square matrix; (0 ≤ n). + a, // Pointer to the n-by-n matrix. + lda, // The leading dimension of `a`. + scratchpad, // Pointer to scratchpad memory to be used by MKL + // routine for storing intermediate results. + scratchpad_size, depends); + } catch (mkl_lapack::exception const &e) { + is_exception_caught = true; + info = e.info(); + if (info < 0) { + error_msg << "Parameter number " << -info + << " had an illegal value."; + } + else if (info == scratchpad_size && e.detail() != 0) { + error_msg + << "Insufficient scratchpad size. Required size is at least " + << e.detail(); + } + else if (info > 0 && e.detail() == 0) { + sycl::free(scratchpad, exec_q); + throw LinAlgError("Matrix is not positive definite."); + } + else { + error_msg << "Unexpected MKL exception caught during getrf() " + "call:\nreason: " + << e.what() << "\ninfo: " << e.info(); + } + } catch (sycl::exception const &e) { + is_exception_caught = true; + error_msg << "Unexpected SYCL exception caught during potrf() call:\n" + << e.what(); + } + + if (is_exception_caught) // an unexpected 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(potrf_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 potrf_event; +} + +std::pair + potrf(sycl::queue q, + dpctl::tensor::usm_ndarray a_array, + const std::int8_t upper_lower, + const std::vector &depends) +{ + const int a_array_nd = a_array.get_ndim(); + + if (a_array_nd != 2) { + throw py::value_error( + "The input array has ndim=" + std::to_string(a_array_nd) + + ", but a 2-dimensional array is expected."); + } + + const py::ssize_t *a_array_shape = a_array.get_shape_raw(); + + if (a_array_shape[0] != a_array_shape[1]) { + throw py::value_error("The input array must be square," + " but got a shape of (" + + std::to_string(a_array_shape[0]) + ", " + + std::to_string(a_array_shape[1]) + ")."); + } + + bool is_a_array_c_contig = a_array.is_c_contiguous(); + if (!is_a_array_c_contig) { + throw py::value_error("The input array " + "must be C-contiguous"); + } + + auto array_types = dpctl_td_ns::usm_ndarray_types(); + int a_array_type_id = + array_types.typenum_to_lookup_id(a_array.get_typenum()); + + potrf_impl_fn_ptr_t potrf_fn = potrf_dispatch_vector[a_array_type_id]; + if (potrf_fn == nullptr) { + throw py::value_error( + "No potrf implementation defined for the provided type " + "of the input matrix."); + } + + char *a_array_data = a_array.get_data(); + const std::int64_t n = a_array_shape[0]; + const std::int64_t lda = std::max(1UL, n); + const oneapi::mkl::uplo uplo_val = + static_cast(upper_lower); + + std::vector host_task_events; + sycl::event potrf_ev = + potrf_fn(q, uplo_val, n, a_array_data, lda, host_task_events, depends); + + sycl::event args_ev = + dpctl::utils::keep_args_alive(q, {a_array}, host_task_events); + + return std::make_pair(args_ev, potrf_ev); +} + +template +struct PotrfContigFactory +{ + fnT get() + { + if constexpr (types::PotrfTypePairSupportFactory::is_defined) { + return potrf_impl; + } + else { + return nullptr; + } + } +}; + +void init_potrf_dispatch_vector(void) +{ + dpctl_td_ns::DispatchVectorBuilder + contig; + contig.populate_dispatch_vector(potrf_dispatch_vector); +} +} // namespace lapack +} // namespace ext +} // namespace backend +} // namespace dpnp diff --git a/dpnp/backend/extensions/lapack/potrf.hpp b/dpnp/backend/extensions/lapack/potrf.hpp new file mode 100644 index 00000000000..f0850b3fd98 --- /dev/null +++ b/dpnp/backend/extensions/lapack/potrf.hpp @@ -0,0 +1,61 @@ +//***************************************************************************** +// Copyright (c) 2024, 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 + potrf(sycl::queue exec_q, + dpctl::tensor::usm_ndarray a_array, + const std::int8_t upper_lower, + const std::vector &depends = {}); + +extern std::pair + potrf_batch(sycl::queue exec_q, + dpctl::tensor::usm_ndarray a_array, + const std::int8_t upper_lower, + const std::int64_t n, + const std::int64_t stride_a, + const std::int64_t batch_size, + const std::vector &depends = {}); + +extern void init_potrf_dispatch_vector(void); +extern void init_potrf_batch_dispatch_vector(void); +} // namespace lapack +} // namespace ext +} // namespace backend +} // namespace dpnp diff --git a/dpnp/backend/extensions/lapack/potrf_batch.cpp b/dpnp/backend/extensions/lapack/potrf_batch.cpp new file mode 100644 index 00000000000..1a36bae4efd --- /dev/null +++ b/dpnp/backend/extensions/lapack/potrf_batch.cpp @@ -0,0 +1,257 @@ +//***************************************************************************** +// Copyright (c) 2024, 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 + +// dpctl tensor headers +#include "utils/memory_overlap.hpp" +#include "utils/type_utils.hpp" + +#include "linalg_exceptions.hpp" +#include "potrf.hpp" +#include "types_matrix.hpp" + +#include "dpnp_utils.hpp" + +namespace dpnp +{ +namespace backend +{ +namespace ext +{ +namespace lapack +{ +namespace mkl_lapack = oneapi::mkl::lapack; +namespace py = pybind11; +namespace type_utils = dpctl::tensor::type_utils; + +typedef sycl::event (*potrf_batch_impl_fn_ptr_t)( + sycl::queue, + const oneapi::mkl::uplo, + const std::int64_t, + char *, + const std::int64_t, + const std::int64_t, + const std::int64_t, + std::vector &, + const std::vector &); + +static potrf_batch_impl_fn_ptr_t + potrf_batch_dispatch_vector[dpctl_td_ns::num_types]; + +template +static sycl::event potrf_batch_impl(sycl::queue exec_q, + const oneapi::mkl::uplo upper_lower, + const std::int64_t n, + char *in_a, + const std::int64_t lda, + const std::int64_t stride_a, + const std::int64_t batch_size, + std::vector &host_task_events, + const std::vector &depends) +{ + type_utils::validate_type_for_device(exec_q); + + T *a = reinterpret_cast(in_a); + + const std::int64_t scratchpad_size = + mkl_lapack::potrf_batch_scratchpad_size(exec_q, upper_lower, n, lda, + stride_a, batch_size); + T *scratchpad = nullptr; + + std::stringstream error_msg; + std::int64_t info = 0; + bool is_exception_caught = false; + + sycl::event potrf_batch_event; + try { + scratchpad = sycl::malloc_device(scratchpad_size, exec_q); + + potrf_batch_event = mkl_lapack::potrf_batch( + exec_q, + upper_lower, // An enumeration value of type oneapi::mkl::uplo: + // oneapi::mkl::uplo::upper for the upper triangular + // part; oneapi::mkl::uplo::lower for the lower + // triangular part. + n, // Order of each square matrix in the batch; (0 ≤ n). + a, // Pointer to the batch of matrices. + lda, // The leading dimension of `a`. + stride_a, // Stride between matrices: Element spacing between + // matrices in `a`. + batch_size, // Total number of matrices in the batch. + scratchpad, // Pointer to scratchpad memory to be used by MKL + // routine for storing intermediate results. + scratchpad_size, depends); + } catch (mkl_lapack::batch_error const &be) { + // Get the indices of matrices within the batch that encountered an + // error + auto error_matrices_ids = be.ids(); + + error_msg + << "Matrix is not positive definite. Errors in matrices with IDs: "; + for (size_t i = 0; i < error_matrices_ids.size(); ++i) { + error_msg << error_matrices_ids[i]; + if (i < error_matrices_ids.size() - 1) { + error_msg << ", "; + } + } + error_msg << "."; + + sycl::free(scratchpad, exec_q); + throw LinAlgError(error_msg.str().c_str()); + } catch (mkl_lapack::exception const &e) { + is_exception_caught = true; + info = e.info(); + + if (info < 0) { + error_msg << "Parameter number " << -info + << " had an illegal value."; + } + else if (info == scratchpad_size && e.detail() != 0) { + error_msg + << "Insufficient scratchpad size. Required size is at least " + << e.detail(); + } + else if (info != 0 && e.detail() == 0) { + error_msg << "Error in batch processing. " + "Number of failed calculations: " + << info; + } + else { + error_msg << "Unexpected MKL exception caught during potrf_batch() " + "call:\nreason: " + << e.what() << "\ninfo: " << e.info(); + } + } catch (sycl::exception const &e) { + is_exception_caught = true; + error_msg + << "Unexpected SYCL exception caught during potrf_batch() call:\n" + << e.what(); + } + + if (is_exception_caught) // an unexpected 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(potrf_batch_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 potrf_batch_event; +} + +std::pair + potrf_batch(sycl::queue q, + dpctl::tensor::usm_ndarray a_array, + const std::int8_t upper_lower, + const std::int64_t n, + const std::int64_t stride_a, + const std::int64_t batch_size, + const std::vector &depends) +{ + const int a_array_nd = a_array.get_ndim(); + + if (a_array_nd < 3) { + throw py::value_error( + "The input array has ndim=" + std::to_string(a_array_nd) + + ", but a 3-dimensional or higher array is expected."); + } + + const py::ssize_t *a_array_shape = a_array.get_shape_raw(); + + if (a_array_shape[a_array_nd - 1] != a_array_shape[a_array_nd - 2]) { + throw py::value_error( + "The last two dimensions of the input array must be square," + " but got a shape of (" + + std::to_string(a_array_shape[a_array_nd - 1]) + ", " + + std::to_string(a_array_shape[a_array_nd - 2]) + ")."); + } + + bool is_a_array_c_contig = a_array.is_c_contiguous(); + if (!is_a_array_c_contig) { + throw py::value_error("The input array " + "must be C-contiguous"); + } + + auto array_types = dpctl_td_ns::usm_ndarray_types(); + int a_array_type_id = + array_types.typenum_to_lookup_id(a_array.get_typenum()); + + potrf_batch_impl_fn_ptr_t potrf_batch_fn = + potrf_batch_dispatch_vector[a_array_type_id]; + if (potrf_batch_fn == nullptr) { + throw py::value_error( + "No potrf_batch implementation defined for the provided type " + "of the input matrix."); + } + + char *a_array_data = a_array.get_data(); + const std::int64_t lda = std::max(1UL, n); + const oneapi::mkl::uplo uplo_val = + static_cast(upper_lower); + + std::vector host_task_events; + sycl::event potrf_batch_ev = + potrf_batch_fn(q, uplo_val, n, a_array_data, lda, stride_a, batch_size, + host_task_events, depends); + + sycl::event args_ev = + dpctl::utils::keep_args_alive(q, {a_array}, host_task_events); + + return std::make_pair(args_ev, potrf_batch_ev); +} + +template +struct PotrfBatchContigFactory +{ + fnT get() + { + if constexpr (types::PotrfBatchTypePairSupportFactory::is_defined) { + return potrf_batch_impl; + } + else { + return nullptr; + } + } +}; + +void init_potrf_batch_dispatch_vector(void) +{ + dpctl_td_ns::DispatchVectorBuilder + contig; + contig.populate_dispatch_vector(potrf_batch_dispatch_vector); +} +} // namespace lapack +} // namespace ext +} // namespace backend +} // namespace dpnp diff --git a/dpnp/backend/extensions/lapack/types_matrix.hpp b/dpnp/backend/extensions/lapack/types_matrix.hpp index 9dca2136f00..0277e630738 100644 --- a/dpnp/backend/extensions/lapack/types_matrix.hpp +++ b/dpnp/backend/extensions/lapack/types_matrix.hpp @@ -142,6 +142,58 @@ struct HeevdTypePairSupportFactory dpctl_td_ns::NotDefinedEntry>::is_defined; }; +/** + * @brief A factory to define pairs of supported types for which + * MKL LAPACK library provides support in oneapi::mkl::lapack::potrf + * function. + * + * @tparam T Type of array containing input matrix, + * as well as the output array for storing the Cholesky factor L. + */ +template +struct PotrfTypePairSupportFactory +{ + static constexpr bool is_defined = std::disjunction< + dpctl_td_ns::TypePairDefinedEntry, + dpctl_td_ns::TypePairDefinedEntry, + dpctl_td_ns::TypePairDefinedEntry, + T, + std::complex>, + dpctl_td_ns::TypePairDefinedEntry, + T, + std::complex>, + // fall-through + dpctl_td_ns::NotDefinedEntry>::is_defined; +}; + +/** + * @brief A factory to define pairs of supported types for which + * MKL LAPACK library provides support in oneapi::mkl::lapack::potrf + * function. + * + * @tparam T Type of array containing input matrices, + * as well as the output arrays for storing the Cholesky factor L. + */ +template +struct PotrfBatchTypePairSupportFactory +{ + static constexpr bool is_defined = std::disjunction< + dpctl_td_ns::TypePairDefinedEntry, + dpctl_td_ns::TypePairDefinedEntry, + dpctl_td_ns::TypePairDefinedEntry, + T, + std::complex>, + dpctl_td_ns::TypePairDefinedEntry, + T, + std::complex>, + // fall-through + dpctl_td_ns::NotDefinedEntry>::is_defined; +}; + /** * @brief A factory to define pairs of supported types for which * MKL LAPACK library provides support in oneapi::mkl::lapack::syevd diff --git a/dpnp/backend/include/dpnp_iface_fptr.hpp b/dpnp/backend/include/dpnp_iface_fptr.hpp index cc0cb9e32a1..e1ee9c76716 100644 --- a/dpnp/backend/include/dpnp_iface_fptr.hpp +++ b/dpnp/backend/include/dpnp_iface_fptr.hpp @@ -88,8 +88,6 @@ enum class DPNPFuncName : size_t DPNP_FN_CBRT, /**< Used in numpy.cbrt() impl */ DPNP_FN_CEIL, /**< Used in numpy.ceil() impl */ DPNP_FN_CHOLESKY, /**< Used in numpy.linalg.cholesky() impl */ - DPNP_FN_CHOLESKY_EXT, /**< Used in numpy.linalg.cholesky() impl, requires - extra parameters */ DPNP_FN_CONJUGATE, /**< Used in numpy.conjugate() impl */ DPNP_FN_CHOOSE, /**< Used in numpy.choose() impl */ DPNP_FN_CHOOSE_EXT, /**< Used in numpy.choose() impl, requires extra diff --git a/dpnp/backend/kernels/dpnp_krnl_linalg.cpp b/dpnp/backend/kernels/dpnp_krnl_linalg.cpp index 8e970884a18..297f963ba28 100644 --- a/dpnp/backend/kernels/dpnp_krnl_linalg.cpp +++ b/dpnp/backend/kernels/dpnp_krnl_linalg.cpp @@ -128,15 +128,6 @@ template void (*dpnp_cholesky_default_c)(void *, void *, const size_t, const size_t) = dpnp_cholesky_c<_DataType>; -template -DPCTLSyclEventRef (*dpnp_cholesky_ext_c)(DPCTLSyclQueueRef, - void *, - void *, - const size_t, - const size_t, - const DPCTLEventVectorRef) = - dpnp_cholesky_c<_DataType>; - template DPCTLSyclEventRef dpnp_det_c(DPCTLSyclQueueRef q_ref, void *array1_in, @@ -860,11 +851,6 @@ void func_map_init_linalg_func(func_map_t &fmap) fmap[DPNPFuncName::DPNP_FN_CHOLESKY][eft_DBL][eft_DBL] = { eft_DBL, (void *)dpnp_cholesky_default_c}; - fmap[DPNPFuncName::DPNP_FN_CHOLESKY_EXT][eft_FLT][eft_FLT] = { - eft_FLT, (void *)dpnp_cholesky_ext_c}; - fmap[DPNPFuncName::DPNP_FN_CHOLESKY_EXT][eft_DBL][eft_DBL] = { - eft_DBL, (void *)dpnp_cholesky_ext_c}; - fmap[DPNPFuncName::DPNP_FN_DET][eft_INT][eft_INT] = { eft_INT, (void *)dpnp_det_default_c}; fmap[DPNPFuncName::DPNP_FN_DET][eft_LNG][eft_LNG] = { diff --git a/dpnp/dpnp_algo/dpnp_algo.pxd b/dpnp/dpnp_algo/dpnp_algo.pxd index 4d59ffcf496..ad36f4137c8 100644 --- a/dpnp/dpnp_algo/dpnp_algo.pxd +++ b/dpnp/dpnp_algo/dpnp_algo.pxd @@ -38,8 +38,6 @@ cdef extern from "dpnp_iface_fptr.hpp" namespace "DPNPFuncName": # need this na DPNP_FN_ARANGE DPNP_FN_ARGSORT DPNP_FN_ARGSORT_EXT - DPNP_FN_CHOLESKY - DPNP_FN_CHOLESKY_EXT DPNP_FN_CHOOSE DPNP_FN_CHOOSE_EXT DPNP_FN_COPY diff --git a/dpnp/linalg/dpnp_algo_linalg.pyx b/dpnp/linalg/dpnp_algo_linalg.pyx index 373f67616fd..6cf3bd3578f 100644 --- a/dpnp/linalg/dpnp_algo_linalg.pyx +++ b/dpnp/linalg/dpnp_algo_linalg.pyx @@ -45,7 +45,6 @@ cimport numpy cimport dpnp.dpnp_utils as utils __all__ = [ - "dpnp_cholesky", "dpnp_cond", "dpnp_eig", "dpnp_eigvals", @@ -67,9 +66,6 @@ ctypedef c_dpctl.DPCTLSyclEventRef(*custom_linalg_1in_1out_func_ptr_t_)(c_dpctl. ctypedef c_dpctl.DPCTLSyclEventRef(*custom_linalg_1in_1out_with_size_func_ptr_t_)(c_dpctl.DPCTLSyclQueueRef, void *, void * , size_t, const c_dpctl.DPCTLEventVectorRef) -ctypedef c_dpctl.DPCTLSyclEventRef(*custom_linalg_1in_1out_with_2size_func_ptr_t_)(c_dpctl.DPCTLSyclQueueRef, - void *, void * , size_t, size_t, - const c_dpctl.DPCTLEventVectorRef) ctypedef c_dpctl.DPCTLSyclEventRef(*custom_linalg_1in_3out_shape_t)(c_dpctl.DPCTLSyclQueueRef, void *, void * , void * , void * , size_t , size_t, const c_dpctl.DPCTLEventVectorRef) @@ -78,43 +74,6 @@ ctypedef c_dpctl.DPCTLSyclEventRef(*custom_linalg_2in_1out_func_ptr_t)(c_dpctl.D const c_dpctl.DPCTLEventVectorRef) -cpdef utils.dpnp_descriptor dpnp_cholesky(utils.dpnp_descriptor input_): - size_ = input_.shape[-1] - - cdef DPNPFuncType param1_type = dpnp_dtype_to_DPNPFuncType(input_.dtype) - - cdef DPNPFuncData kernel_data = get_dpnp_function_ptr(DPNP_FN_CHOLESKY_EXT, param1_type, param1_type) - - input_obj = input_.get_array() - - # create result array with type given by FPTR data - cdef utils.dpnp_descriptor result = utils.create_output_descriptor(input_.shape, - kernel_data.return_type, - None, - device=input_obj.sycl_device, - usm_type=input_obj.usm_type, - sycl_queue=input_obj.sycl_queue) - - result_sycl_queue = result.get_array().sycl_queue - - cdef c_dpctl.SyclQueue q = result_sycl_queue - cdef c_dpctl.DPCTLSyclQueueRef q_ref = q.get_queue_ref() - - cdef custom_linalg_1in_1out_with_2size_func_ptr_t_ func = kernel_data.ptr - - cdef c_dpctl.DPCTLSyclEventRef event_ref = func(q_ref, - input_.get_data(), - result.get_data(), - input_.size, - size_, - NULL) # dep_events_ref - - with nogil: c_dpctl.DPCTLEvent_WaitAndThrow(event_ref) - c_dpctl.DPCTLEvent_Delete(event_ref) - - return result - - cpdef object dpnp_cond(object input, object p): if p in ('f', 'fro'): # TODO: change order='K' when support is implemented diff --git a/dpnp/linalg/dpnp_iface_linalg.py b/dpnp/linalg/dpnp_iface_linalg.py index 008a7a2212b..2f7420d1329 100644 --- a/dpnp/linalg/dpnp_iface_linalg.py +++ b/dpnp/linalg/dpnp_iface_linalg.py @@ -47,6 +47,7 @@ from .dpnp_utils_linalg import ( check_stacked_2d, check_stacked_square, + dpnp_cholesky, dpnp_det, dpnp_eigh, dpnp_slogdet, @@ -72,52 +73,64 @@ ] -def cholesky(input): +def cholesky(a, upper=False): """ Cholesky decomposition. - Return the Cholesky decomposition, `L * L.H`, of the square matrix `input`, - where `L` is lower-triangular and .H is the conjugate transpose operator - (which is the ordinary transpose if `input` is real-valued). `input` must be - Hermitian (symmetric if real-valued) and positive-definite. No - checking is performed to verify whether `a` is Hermitian or not. - In addition, only the lower-triangular and diagonal elements of `input` - are used. Only `L` is actually returned. + Return the lower or upper Cholesky decomposition, ``L * L.H`` or + ``U.H * U``, of the square matrix ``a``, where ``L`` is lower-triangular, + ``U`` is upper-triangular, and ``.H`` is the conjugate transpose operator + (which is the ordinary transpose if ``a`` is real-valued). ``a`` must be + Hermitian (symmetric if real-valued) and positive-definite. No checking is + performed to verify whether ``a`` is Hermitian or not. In addition, only + the lower or upper-triangular and diagonal elements of ``a`` are used. + Only ``L`` or ``U`` is actually returned. + + For full documentation refer to :obj:`numpy.linalg.cholesky`. Parameters ---------- - input : (..., M, M) array_like + a : (..., M, M) {dpnp.ndarray, usm_ndarray} Hermitian (symmetric if all elements are real), positive-definite input matrix. + upper : bool, optional + If ``True``, the result must be the upper-triangular Cholesky factor. + If ``False``, the result must be the lower-triangular Cholesky factor. + Default: ``False``. Returns ------- - L : (..., M, M) array_like - Upper or lower-triangular Cholesky factor of `input`. Returns a - matrix object if `input` is a matrix object. + L : (..., M, M) dpnp.ndarray + Lower or upper-triangular Cholesky factor of `a`. + + Examples + -------- + >>> import dpnp as np + >>> A = np.array([[1.0, 2.0],[2.0, 5.0]]) + >>> A + array([[1., 2.], + [2., 5.]]) + >>> L = np.linalg.cholesky(A) + >>> L + array([[1., 0.], + [2., 1.]]) + >>> np.dot(L, L.T.conj()) # verify that L * L.H = A + array([[1., 2.], + [2., 5.]]) + + The upper-triangular Cholesky factor can also be obtained: + + >>> np.linalg.cholesky(A, upper=True) + array([[ 1.+0.j, -0.-2.j], + [ 0.+0.j, 1.+0.j]] + """ - x1_desc = dpnp.get_dpnp_descriptor(input, copy_when_nondefault_queue=False) - if x1_desc: - if x1_desc.shape[-1] != x1_desc.shape[-2]: - pass - else: - if input.dtype == dpnp.int32 or input.dtype == dpnp.int64: - dev = x1_desc.get_array().sycl_device - if dev.has_aspect_fp64: - dtype = dpnp.float64 - else: - dtype = dpnp.float32 - # TODO memory copy. needs to move into DPNPC - input_ = dpnp.get_dpnp_descriptor( - dpnp.astype(input, dtype=dtype), - copy_when_nondefault_queue=False, - ) - else: - input_ = x1_desc - return dpnp_cholesky(input_).get_pyobj() - - return call_origin(numpy.linalg.cholesky, input) + dpnp.check_supported_arrays_type(a) + check_stacked_2d(a) + check_stacked_square(a) + + return dpnp_cholesky(a, upper=upper) def cond(input, p=None): diff --git a/dpnp/linalg/dpnp_utils_linalg.py b/dpnp/linalg/dpnp_utils_linalg.py index 62d20cb51a8..40159ac02bc 100644 --- a/dpnp/linalg/dpnp_utils_linalg.py +++ b/dpnp/linalg/dpnp_utils_linalg.py @@ -34,6 +34,7 @@ __all__ = [ "check_stacked_2d", "check_stacked_square", + "dpnp_cholesky", "dpnp_det", "dpnp_eigh", "dpnp_slogdet", @@ -440,6 +441,130 @@ def _lu_factor(a, res_type): return (a_h, ipiv_h, dev_info_array) +def dpnp_cholesky_batch(a, upper_lower, res_type): + """ + dpnp_cholesky_batch(a, upper_lower, res_type) + + Return the batched Cholesky decomposition of `a` array. + + """ + + a_sycl_queue = a.sycl_queue + a_usm_type = a.usm_type + + n = a.shape[-2] + + orig_shape = a.shape + # get 3d input arrays by reshape + a = a.reshape(-1, n, n) + batch_size = a.shape[0] + a_usm_arr = dpnp.get_usm_ndarray(a) + + # `a` must be copied because potrf_batch destroys the input matrix + a_h = dpnp.empty_like(a, order="C", dtype=res_type, usm_type=a_usm_type) + + # use DPCTL tensor function to fill the сopy of the input array + # from the input array + a_ht_copy_ev, a_copy_ev = ti._copy_usm_ndarray_into_usm_ndarray( + src=a_usm_arr, dst=a_h.get_array(), sycl_queue=a_sycl_queue + ) + + a_stride = a_h.strides[0] + + # Call the LAPACK extension function _potrf_batch + # to computes the Cholesky decomposition of a batch of + # symmetric positive-definite matrices + ht_lapack_ev, _ = li._potrf_batch( + a_sycl_queue, + a_h.get_array(), + upper_lower, + n, + a_stride, + batch_size, + [a_copy_ev], + ) + + ht_lapack_ev.wait() + a_ht_copy_ev.wait() + + # Get upper or lower-triangular matrix part as per `upper_lower` value + # upper_lower is 0 (lower) or 1 (upper) + if upper_lower: + a_h = dpnp.triu(a_h.reshape(orig_shape)) + else: + a_h = dpnp.tril(a_h.reshape(orig_shape)) + + return a_h + + +def dpnp_cholesky(a, upper): + """ + dpnp_cholesky(a, upper) + + Return the Cholesky decomposition of `a` array. + + """ + + a_sycl_queue = a.sycl_queue + a_usm_type = a.usm_type + + res_type = _common_type(a) + + a_shape = a.shape + + if a.size == 0: + return dpnp.empty( + a_shape, + dtype=res_type, + usm_type=a_usm_type, + sycl_queue=a_sycl_queue, + ) + + # Set `uplo` value for `potrf` and `potrf_batch` function based on the boolean input `upper`. + # In oneMKL, `uplo` value of 1 is equivalent to oneapi::mkl::uplo::lower + # and `uplo` value of 0 is equivalent to oneapi::mkl::uplo::upper. + # However, we adjust this logic based on the array's memory layout. + # Note: lower for row-major (which is used here) is upper for column-major layout. + # Reference: comment from tbmkl/tests/lapack/unit/dpcpp/potrf_usm/potrf_usm.cpp + # This means that if `upper` is False (lower triangular), + # we actually use oneapi::mkl::uplo::upper (0) for the row-major layout, and vice versa. + upper_lower = int(upper) + + if a.ndim > 2: + return dpnp_cholesky_batch(a, upper_lower, res_type) + + a_usm_arr = dpnp.get_usm_ndarray(a) + + # `a` must be copied because potrf destroys the input matrix + a_h = dpnp.empty_like(a, order="C", dtype=res_type, usm_type=a_usm_type) + + # use DPCTL tensor function to fill the сopy of the input array + # from the input array + a_ht_copy_ev, a_copy_ev = ti._copy_usm_ndarray_into_usm_ndarray( + src=a_usm_arr, dst=a_h.get_array(), sycl_queue=a_sycl_queue + ) + + # Call the LAPACK extension function _potrf + # to computes the Cholesky decomposition + ht_lapack_ev, _ = li._potrf( + a_sycl_queue, + a_h.get_array(), + upper_lower, + [a_copy_ev], + ) + + ht_lapack_ev.wait() + a_ht_copy_ev.wait() + + # Get upper or lower-triangular matrix part as per `upper` value + if upper: + a_h = dpnp.triu(a_h) + else: + a_h = dpnp.tril(a_h) + + return a_h + + def dpnp_det(a): """ dpnp_det(a) diff --git a/tests/test_linalg.py b/tests/test_linalg.py index cb780abd925..5f38421c6ec 100644 --- a/tests/test_linalg.py +++ b/tests/test_linalg.py @@ -4,6 +4,7 @@ from numpy.testing import assert_allclose, assert_array_equal, assert_raises import dpnp as inp +from tests.third_party.cupy import testing from .helper import ( assert_dtype_allclose, @@ -44,44 +45,159 @@ def vvsort(val, vec, size, xp): vec[:, imax] = temp -@pytest.mark.parametrize( - "array", - [ - [[[1, -2], [2, 5]]], - [[[1.0, -2.0], [2.0, 5.0]]], - [[[1.0, -2.0], [2.0, 5.0]], [[1.0, -2.0], [2.0, 5.0]]], - ], - ids=[ - "[[[1, -2], [2, 5]]]", - "[[[1., -2.], [2., 5.]]]", - "[[[1., -2.], [2., 5.]], [[1., -2.], [2., 5.]]]", - ], -) -def test_cholesky(array): - a = numpy.array(array) - ia = inp.array(a) - result = inp.linalg.cholesky(ia) - expected = numpy.linalg.cholesky(a) - assert_array_equal(expected, result) +class TestCholesky: + @pytest.mark.parametrize( + "array", + [ + [[1, 2], [2, 5]], + [[[5, 2], [2, 6]], [[7, 3], [3, 8]], [[3, 1], [1, 4]]], + [ + [[[5, 2], [2, 5]], [[6, 3], [3, 6]]], + [[[7, 2], [2, 7]], [[8, 3], [3, 8]]], + ], + ], + ids=[ + "2D_array", + "3D_array", + "4D_array", + ], + ) + @pytest.mark.parametrize("dtype", get_all_dtypes(no_bool=True)) + def test_cholesky(self, array, dtype): + a = numpy.array(array, dtype=dtype) + ia = inp.array(a) + result = inp.linalg.cholesky(ia) + expected = numpy.linalg.cholesky(a) + assert_dtype_allclose(result, expected) + + @pytest.mark.parametrize( + "array", + [ + [[1, 2], [2, 5]], + [[[5, 2], [2, 6]], [[7, 3], [3, 8]], [[3, 1], [1, 4]]], + [ + [[[5, 2], [2, 5]], [[6, 3], [3, 6]]], + [[[7, 2], [2, 7]], [[8, 3], [3, 8]]], + ], + ], + ids=[ + "2D_array", + "3D_array", + "4D_array", + ], + ) + @pytest.mark.parametrize("dtype", get_all_dtypes(no_bool=True)) + def test_cholesky_upper(self, array, dtype): + ia = inp.array(array, dtype=dtype) + result = inp.linalg.cholesky(ia, upper=True) + + if ia.ndim > 2: + n = ia.shape[-1] + ia_reshaped = ia.reshape(-1, n, n) + res_reshaped = result.reshape(-1, n, n) + batch_size = ia_reshaped.shape[0] + for idx in range(batch_size): + # Reconstruct the matrix using the Cholesky decomposition result + if inp.issubdtype(dtype, inp.complexfloating): + reconstructed = ( + res_reshaped[idx].T.conj() @ res_reshaped[idx] + ) + else: + reconstructed = res_reshaped[idx].T @ res_reshaped[idx] + assert_dtype_allclose( + reconstructed, ia_reshaped[idx], check_type=False + ) + else: + # Reconstruct the matrix using the Cholesky decomposition result + if inp.issubdtype(dtype, inp.complexfloating): + reconstructed = result.T.conj() @ result + else: + reconstructed = result.T @ result + assert_dtype_allclose(reconstructed, ia, check_type=False) + + # upper parameter support will be added in numpy 2.0 version + @testing.with_requires("numpy>=2.0") + @pytest.mark.parametrize( + "array", + [ + [[1, 2], [2, 5]], + [[[5, 2], [2, 6]], [[7, 3], [3, 8]], [[3, 1], [1, 4]]], + [ + [[[5, 2], [2, 5]], [[6, 3], [3, 6]]], + [[[7, 2], [2, 7]], [[8, 3], [3, 8]]], + ], + ], + ids=[ + "2D_array", + "3D_array", + "4D_array", + ], + ) + @pytest.mark.parametrize("dtype", get_all_dtypes(no_bool=True)) + def test_cholesky_upper_numpy(self, array, dtype): + a = numpy.array(array, dtype=dtype) + ia = inp.array(a) + result = inp.linalg.cholesky(ia, upper=True) + expected = numpy.linalg.cholesky(a, upper=True) + assert_dtype_allclose(result, expected) + def test_cholesky_strides(self): + a_np = numpy.array( + [ + [5, 2, 0, 0, 1], + [2, 6, 0, 0, 2], + [0, 0, 7, 0, 0], + [0, 0, 0, 4, 0], + [1, 2, 0, 0, 5], + ] + ) -@pytest.mark.parametrize( - "shape", - [ - (0, 0), - (3, 0, 0), - ], - ids=[ - "(0, 0)", - "(3, 0, 0)", - ], -) -def test_cholesky_0D(shape): - a = numpy.empty(shape) - ia = inp.array(a) - result = inp.linalg.cholesky(ia) - expected = numpy.linalg.cholesky(a) - assert_array_equal(expected, result) + a_dp = inp.array(a_np) + + # positive strides + expected = numpy.linalg.cholesky(a_np[::2, ::2]) + result = inp.linalg.cholesky(a_dp[::2, ::2]) + assert_allclose(expected, result, rtol=1e-3, atol=1e-4) + + # negative strides + expected = numpy.linalg.cholesky(a_np[::-2, ::-2]) + result = inp.linalg.cholesky(a_dp[::-2, ::-2]) + assert_allclose(expected, result, rtol=1e-3, atol=1e-4) + + @pytest.mark.parametrize( + "shape", + [ + (0, 0), + (3, 0, 0), + (0, 2, 2), + ], + ids=[ + "(0, 0)", + "(3, 0, 0)", + "(0, 2, 2)", + ], + ) + def test_cholesky_empty(self, shape): + a = numpy.empty(shape) + ia = inp.array(a) + result = inp.linalg.cholesky(ia) + expected = numpy.linalg.cholesky(a) + assert_array_equal(expected, result) + + def test_cholesky_errors(self): + a_dp = inp.array([[1, 2], [2, 5]], dtype="float32") + + # unsupported type + a_np = inp.asnumpy(a_dp) + assert_raises(TypeError, inp.linalg.cholesky, a_np) + + # a.ndim < 2 + a_dp_ndim_1 = a_dp.flatten() + assert_raises(inp.linalg.LinAlgError, inp.linalg.cholesky, a_dp_ndim_1) + + # a is not square + a_dp = inp.ones((2, 3)) + assert_raises(inp.linalg.LinAlgError, inp.linalg.cholesky, a_dp) @pytest.mark.parametrize( diff --git a/tests/test_random_state.py b/tests/test_random_state.py index 4771eadc42e..70940501d2e 100644 --- a/tests/test_random_state.py +++ b/tests/test_random_state.py @@ -491,7 +491,7 @@ def test_rng_zero_and_extremes(self): sycl_device = dpctl.SyclQueue().sycl_device if sycl_device.has_aspect_gpu and not sycl_device.has_aspect_fp64: - # TODO: discuss with opneMKL + # TODO: discuss with oneMKL pytest.skip( f"Due to some reason, oneMKL wrongly returns high value instead of low" ) diff --git a/tests/test_sycl_queue.py b/tests/test_sycl_queue.py index fd09d08b5be..0e66a3214bd 100644 --- a/tests/test_sycl_queue.py +++ b/tests/test_sycl_queue.py @@ -960,19 +960,31 @@ def test_fft_rfft(type, shape, device): assert_sycl_queue_equal(result_queue, expected_queue) +@pytest.mark.parametrize( + "data, is_empty", + [ + ([[1, -2], [2, 5]], False), + ([[[1, -2], [2, 5]], [[1, -2], [2, 5]]], False), + ((0, 0), True), + ((3, 0, 0), True), + ], + ids=["2D", "3D", "Empty_2D", "Empty_3D"], +) @pytest.mark.parametrize( "device", valid_devices, ids=[device.filter_string for device in valid_devices], ) -def test_cholesky(device): - data = [[[1.0, -2.0], [2.0, 5.0]], [[1.0, -2.0], [2.0, 5.0]]] - numpy_data = numpy.array(data) - dpnp_data = dpnp.array(data, device=device) +def test_cholesky(data, is_empty, device): + if is_empty: + numpy_data = numpy.empty(data, dtype=dpnp.default_float_type(device)) + else: + numpy_data = numpy.array(data, dtype=dpnp.default_float_type(device)) + dpnp_data = dpnp.array(numpy_data, device=device) result = dpnp.linalg.cholesky(dpnp_data) expected = numpy.linalg.cholesky(numpy_data) - assert_array_equal(expected, result) + assert_dtype_allclose(result, expected) expected_queue = dpnp_data.get_array().sycl_queue result_queue = result.get_array().sycl_queue diff --git a/tests/test_usm_type.py b/tests/test_usm_type.py index fa4a14090b3..d2bf335effb 100644 --- a/tests/test_usm_type.py +++ b/tests/test_usm_type.py @@ -558,6 +558,28 @@ def test_take(func, usm_type_x, usm_type_ind): assert z.usm_type == du.get_coerced_usm_type([usm_type_x, usm_type_ind]) +@pytest.mark.parametrize( + "data, is_empty", + [ + ([[1, -2], [2, 5]], False), + ([[[1, -2], [2, 5]], [[1, -2], [2, 5]]], False), + ((0, 0), True), + ((3, 0, 0), True), + ], + ids=["2D", "3D", "Empty_2D", "Empty_3D"], +) +@pytest.mark.parametrize("usm_type", list_of_usm_types, ids=list_of_usm_types) +def test_cholesky(data, is_empty, usm_type): + if is_empty: + x = dp.empty(data, dtype=dp.default_float_type(), usm_type=usm_type) + else: + x = dp.array(data, dtype=dp.default_float_type(), usm_type=usm_type) + + result = dp.linalg.cholesky(x) + + assert x.usm_type == result.usm_type + + @pytest.mark.parametrize("usm_type", list_of_usm_types, ids=list_of_usm_types) def test_indices(usm_type): x = dp.indices((2,), usm_type=usm_type) diff --git a/tests/third_party/cupy/linalg_tests/test_decomposition.py b/tests/third_party/cupy/linalg_tests/test_decomposition.py new file mode 100644 index 00000000000..42bcf122ff4 --- /dev/null +++ b/tests/third_party/cupy/linalg_tests/test_decomposition.py @@ -0,0 +1,137 @@ +import unittest + +import numpy +import pytest + +import dpnp as cupy +from tests.helper import has_support_aspect64, is_cpu_device +from tests.third_party.cupy import testing + + +def random_matrix(shape, dtype, scale, sym=False): + m, n = shape[-2:] + dtype = numpy.dtype(dtype) + assert dtype.kind in "iufc" + low_s, high_s = scale + bias = None + if dtype.kind in "iu": + # For an m \times n matrix M whose element is in [-0.5, 0.5], it holds + # (singular value of M) <= \sqrt{mn} / 2 + err = numpy.sqrt(m * n) / 2.0 + low_s += err + high_s -= err + if dtype.kind in "u": + assert sym, ( + "generating nonsymmetric matrix with uint cells is not" + " supported." + ) + # (singular value of numpy.ones((m, n))) <= \sqrt{mn} + high_s = bias = high_s / (1 + numpy.sqrt(m * n)) + assert low_s <= high_s + a = numpy.random.standard_normal(shape) + if dtype.kind == "c": + a = a + 1j * numpy.random.standard_normal(shape) + u, s, vh = numpy.linalg.svd(a) + if sym: + assert m == n + vh = u.conj().swapaxes(-1, -2) + new_s = numpy.random.uniform(low_s, high_s, s.shape) + new_a = numpy.einsum("...ij,...j,...jk->...ik", u, new_s, vh) + if bias is not None: + new_a += bias + if dtype.kind in "iu": + new_a = numpy.rint(new_a) + return new_a.astype(dtype) + + +class TestCholeskyDecomposition: + @testing.numpy_cupy_allclose(atol=1e-3, type_check=has_support_aspect64()) + def check_L(self, array, xp): + a = xp.asarray(array) + return xp.linalg.cholesky(a) + + @testing.for_dtypes( + [ + numpy.int32, + numpy.int64, + numpy.uint32, + numpy.uint64, + numpy.float32, + numpy.float64, + numpy.complex64, + numpy.complex128, + ] + ) + def test_decomposition(self, dtype): + # A positive definite matrix + A = random_matrix((5, 5), dtype, scale=(10, 10000), sym=True) + self.check_L(A) + # np.linalg.cholesky only uses a lower triangle of an array + self.check_L(numpy.array([[1, 2], [1, 9]], dtype)) + + @testing.for_dtypes( + [ + numpy.int32, + numpy.int64, + numpy.uint32, + numpy.uint64, + numpy.float32, + numpy.float64, + numpy.complex64, + numpy.complex128, + ] + ) + def test_batched_decomposition(self, dtype): + Ab1 = random_matrix((3, 5, 5), dtype, scale=(10, 10000), sym=True) + self.check_L(Ab1) + Ab2 = random_matrix((2, 2, 5, 5), dtype, scale=(10, 10000), sym=True) + self.check_L(Ab2) + + @pytest.mark.parametrize( + "shape", + [ + # empty square + (0, 0), + (3, 0, 0), + # empty batch + (2, 0, 3, 4, 4), + ], + ) + @testing.for_dtypes( + [ + numpy.int32, + numpy.uint16, + numpy.float32, + numpy.float64, + numpy.complex64, + numpy.complex128, + ] + ) + @testing.numpy_cupy_allclose(type_check=has_support_aspect64()) + def test_empty(self, shape, xp, dtype): + a = xp.empty(shape, dtype=dtype) + return xp.linalg.cholesky(a) + + +class TestCholeskyInvalid(unittest.TestCase): + def check_L(self, array): + for xp in (numpy, cupy): + a = xp.asarray(array) + with pytest.raises(xp.linalg.LinAlgError): + xp.linalg.cholesky(a) + + # TODO: remove skipif when MKLD-16626 is resolved + @pytest.mark.skipif(is_cpu_device(), reason="MKLD-16626") + @testing.for_dtypes( + [ + numpy.int32, + numpy.int64, + numpy.uint32, + numpy.uint64, + numpy.float32, + numpy.float64, + ] + ) + def test_decomposition(self, dtype): + A = numpy.array([[1, -2], [-2, 1]]).astype(dtype) + self.check_L(A)