diff --git a/.github/workflows/conda-package.yml b/.github/workflows/conda-package.yml index 9b6b093801ac..b633b4bdc880 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_norms.py third_party/cupy/linalg_tests/test_product.py third_party/cupy/linalg_tests/test_solve.py third_party/cupy/logic_tests/test_comparison.py diff --git a/dpnp/backend/extensions/lapack/CMakeLists.txt b/dpnp/backend/extensions/lapack/CMakeLists.txt index d224c623c8cb..77ba9cd84ea0 100644 --- a/dpnp/backend/extensions/lapack/CMakeLists.txt +++ b/dpnp/backend/extensions/lapack/CMakeLists.txt @@ -28,6 +28,8 @@ set(python_module_name _lapack_impl) set(_module_src ${CMAKE_CURRENT_SOURCE_DIR}/lapack_py.cpp ${CMAKE_CURRENT_SOURCE_DIR}/gesv.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/getrf.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/getrf_batch.cpp ${CMAKE_CURRENT_SOURCE_DIR}/heevd.cpp ${CMAKE_CURRENT_SOURCE_DIR}/syevd.cpp ) diff --git a/dpnp/backend/extensions/lapack/gesv.cpp b/dpnp/backend/extensions/lapack/gesv.cpp index 72e5aa806714..781356a2f52a 100644 --- a/dpnp/backend/extensions/lapack/gesv.cpp +++ b/dpnp/backend/extensions/lapack/gesv.cpp @@ -84,7 +84,7 @@ static sycl::event gesv_impl(sycl::queue exec_q, std::stringstream error_msg; std::int64_t info = 0; - bool sycl_exception_caught = false; + bool is_exception_caught = false; sycl::event gesv_event; try { @@ -106,12 +106,18 @@ static sycl::event gesv_impl(sycl::queue exec_q, // 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) { T host_U; exec_q.memcpy(&host_U, &a[(info - 1) * lda + info - 1], sizeof(T)) @@ -131,23 +137,18 @@ static sycl::event gesv_impl(sycl::queue exec_q, << e.what() << "\ninfo: " << e.info(); } } - else if (info == scratchpad_size && e.detail() != 0) { - error_msg - << "Insufficient scratchpad size. Required size is at least " - << e.detail(); - } else { error_msg << "Unexpected MKL exception caught during gesv() " "call:\nreason: " << e.what() << "\ninfo: " << e.info(); } } catch (sycl::exception const &e) { + is_exception_caught = true; error_msg << "Unexpected SYCL exception caught during gesv() call:\n" << e.what(); - sycl_exception_caught = true; } - if (info != 0 || sycl_exception_caught) // an unexpected error occurs + if (is_exception_caught) // an unexpected error occurs { if (scratchpad != nullptr) { sycl::free(scratchpad, exec_q); diff --git a/dpnp/backend/extensions/lapack/getrf.cpp b/dpnp/backend/extensions/lapack/getrf.cpp new file mode 100644 index 000000000000..faebb9f7d42f --- /dev/null +++ b/dpnp/backend/extensions/lapack/getrf.cpp @@ -0,0 +1,256 @@ +//***************************************************************************** +// 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 + +// dpctl tensor headers +#include "utils/memory_overlap.hpp" +#include "utils/type_utils.hpp" + +#include "getrf.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 (*getrf_impl_fn_ptr_t)(sycl::queue, + const std::int64_t, + char *, + std::int64_t, + std::int64_t *, + py::list, + std::vector &, + const std::vector &); + +static getrf_impl_fn_ptr_t getrf_dispatch_vector[dpctl_td_ns::num_types]; + +template +static sycl::event getrf_impl(sycl::queue exec_q, + const std::int64_t n, + char *in_a, + std::int64_t lda, + std::int64_t *ipiv, + py::list dev_info, + 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::getrf_scratchpad_size(exec_q, n, n, lda); + T *scratchpad = nullptr; + + std::stringstream error_msg; + std::int64_t info = 0; + bool is_exception_caught = false; + + sycl::event getrf_event; + try { + scratchpad = sycl::malloc_device(scratchpad_size, exec_q); + + getrf_event = mkl_lapack::getrf( + exec_q, + n, // The order of the square matrix A (0 ≤ n). + // It must be a non-negative integer. + n, // The number of columns in the square matrix A (0 ≤ n). + // It must be a non-negative integer. + a, // Pointer to the square matrix A (n x n). + lda, // The leading dimension of matrix A. + // It must be at least max(1, n). + ipiv, // Pointer to the output array of pivot indices. + 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) { + // Store the positive 'info' value in the first element of + // 'dev_info'. This indicates that the factorization has been + // completed, but the factor U (upper triangular matrix) is exactly + // singular. The 'info' value here is the index of the first zero + // element in the diagonal of U. + is_exception_caught = false; + dev_info[0] = info; + } + 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 getrf() 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(getrf_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 getrf_event; +} + +std::pair + getrf(sycl::queue exec_q, + dpctl::tensor::usm_ndarray a_array, + dpctl::tensor::usm_ndarray ipiv_array, + py::list dev_info, + const std::vector &depends) +{ + const int a_array_nd = a_array.get_ndim(); + const int ipiv_array_nd = ipiv_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."); + } + + if (ipiv_array_nd != 1) { + throw py::value_error("The array of pivot indices has ndim=" + + std::to_string(ipiv_array_nd) + + ", but a 1-dimensional array is expected."); + } + + // check compatibility of execution queue and allocation queue + if (!dpctl::utils::queues_are_compatible(exec_q, {a_array, ipiv_array})) { + throw py::value_error( + "Execution queue is not compatible with allocation queues"); + } + + auto const &overlap = dpctl::tensor::overlap::MemoryOverlap(); + if (overlap(a_array, ipiv_array)) { + throw py::value_error("The input array and the array of pivot indices " + "are overlapping segments of memory"); + } + + bool is_a_array_c_contig = a_array.is_c_contiguous(); + bool is_ipiv_array_c_contig = ipiv_array.is_c_contiguous(); + if (!is_a_array_c_contig) { + throw py::value_error("The input array " + "must be C-contiguous"); + } + if (!is_ipiv_array_c_contig) { + throw py::value_error("The array of pivot indices " + "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()); + + getrf_impl_fn_ptr_t getrf_fn = getrf_dispatch_vector[a_array_type_id]; + if (getrf_fn == nullptr) { + throw py::value_error( + "No getrf implementation defined for the provided type " + "of the input matrix."); + } + + auto ipiv_types = dpctl_td_ns::usm_ndarray_types(); + int ipiv_array_type_id = + ipiv_types.typenum_to_lookup_id(ipiv_array.get_typenum()); + + if (ipiv_array_type_id != static_cast(dpctl_td_ns::typenum_t::INT64)) { + throw py::value_error("The type of 'ipiv_array' must be int64."); + } + + const std::int64_t n = a_array.get_shape_raw()[0]; + + char *a_array_data = a_array.get_data(); + const std::int64_t lda = std::max(1UL, n); + + char *ipiv_array_data = ipiv_array.get_data(); + std::int64_t *d_ipiv = reinterpret_cast(ipiv_array_data); + + std::vector host_task_events; + sycl::event getrf_ev = getrf_fn(exec_q, n, a_array_data, lda, d_ipiv, + dev_info, host_task_events, depends); + + sycl::event args_ev = dpctl::utils::keep_args_alive( + exec_q, {a_array, ipiv_array}, host_task_events); + + return std::make_pair(args_ev, getrf_ev); +} + +template +struct GetrfContigFactory +{ + fnT get() + { + if constexpr (types::GetrfTypePairSupportFactory::is_defined) { + return getrf_impl; + } + else { + return nullptr; + } + } +}; + +void init_getrf_dispatch_vector(void) +{ + dpctl_td_ns::DispatchVectorBuilder + contig; + contig.populate_dispatch_vector(getrf_dispatch_vector); +} +} // namespace lapack +} // namespace ext +} // namespace backend +} // namespace dpnp diff --git a/dpnp/backend/extensions/lapack/getrf.hpp b/dpnp/backend/extensions/lapack/getrf.hpp new file mode 100644 index 000000000000..d8ea425f6804 --- /dev/null +++ b/dpnp/backend/extensions/lapack/getrf.hpp @@ -0,0 +1,64 @@ +//***************************************************************************** +// 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 + getrf(sycl::queue exec_q, + dpctl::tensor::usm_ndarray a_array, + dpctl::tensor::usm_ndarray ipiv_array, + py::list dev_info, + const std::vector &depends = {}); + +extern std::pair + getrf_batch(sycl::queue exec_q, + dpctl::tensor::usm_ndarray a_array, + dpctl::tensor::usm_ndarray ipiv_array, + py::list dev_info, + std::int64_t n, + std::int64_t stride_a, + std::int64_t stride_ipiv, + std::int64_t batch_size, + const std::vector &depends = {}); + +extern void init_getrf_dispatch_vector(void); +extern void init_getrf_batch_dispatch_vector(void); +} // namespace lapack +} // namespace ext +} // namespace backend +} // namespace dpnp diff --git a/dpnp/backend/extensions/lapack/getrf_batch.cpp b/dpnp/backend/extensions/lapack/getrf_batch.cpp new file mode 100644 index 000000000000..0944d7437aed --- /dev/null +++ b/dpnp/backend/extensions/lapack/getrf_batch.cpp @@ -0,0 +1,278 @@ +//***************************************************************************** +// 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 + +// dpctl tensor headers +#include "utils/memory_overlap.hpp" +#include "utils/type_utils.hpp" + +#include "getrf.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 (*getrf_batch_impl_fn_ptr_t)( + sycl::queue, + std::int64_t, + char *, + std::int64_t, + std::int64_t, + std::int64_t *, + std::int64_t, + std::int64_t, + py::list, + std::vector &, + const std::vector &); + +static getrf_batch_impl_fn_ptr_t + getrf_batch_dispatch_vector[dpctl_td_ns::num_types]; + +template +static sycl::event getrf_batch_impl(sycl::queue exec_q, + std::int64_t n, + char *in_a, + std::int64_t lda, + std::int64_t stride_a, + std::int64_t *ipiv, + std::int64_t stride_ipiv, + std::int64_t batch_size, + py::list dev_info, + 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::getrf_batch_scratchpad_size(exec_q, n, n, lda, stride_a, + stride_ipiv, batch_size); + T *scratchpad = nullptr; + + std::stringstream error_msg; + std::int64_t info = 0; + bool is_exception_caught = false; + + sycl::event getrf_batch_event; + try { + scratchpad = sycl::malloc_device(scratchpad_size, exec_q); + + getrf_batch_event = mkl_lapack::getrf_batch( + exec_q, + n, // The order of each square matrix in the batch; (0 ≤ n). + // It must be a non-negative integer. + n, // The number of columns in each matrix in the batch; (0 ≤ n). + // It must be a non-negative integer. + a, // Pointer to the batch of square matrices, each of size (n x n). + lda, // The leading dimension of each matrix in the batch. + stride_a, // Stride between consecutive matrices in the batch. + ipiv, // Pointer to the array of pivot indices for each matrix in + // the batch. + stride_ipiv, // Stride between pivot indices: Spacing between pivot + // arrays in 'ipiv'. + batch_size, // Stride between pivot index arrays 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(); + // Get the indices of the first zero diagonal elements of these matrices + auto error_info = be.exceptions(); + + for (size_t i = 0; i < error_matrices_ids.size(); ++i) { + // Assign the index of the first zero diagonal element in each + // error matrix to the corresponding index in 'dev_info' + dev_info[error_matrices_ids[i]] = error_info[i]; + } + } 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 { + error_msg << "Unexpected MKL exception caught during getrf_batch() " + "call:\nreason: " + << e.what() << "\ninfo: " << e.info(); + } + } catch (sycl::exception const &e) { + is_exception_caught = true; + error_msg + << "Unexpected SYCL exception caught during getrf_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(getrf_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 getrf_batch_event; +} + +std::pair + getrf_batch(sycl::queue exec_q, + dpctl::tensor::usm_ndarray a_array, + dpctl::tensor::usm_ndarray ipiv_array, + py::list dev_info, + std::int64_t n, + std::int64_t stride_a, + std::int64_t stride_ipiv, + std::int64_t batch_size, + const std::vector &depends) +{ + const int a_array_nd = a_array.get_ndim(); + const int ipiv_array_nd = ipiv_array.get_ndim(); + + if (a_array_nd < 3) { + throw py::value_error( + "The input array has ndim=" + std::to_string(a_array_nd) + + ", but an array with ndim >= 3 is expected."); + } + + if (ipiv_array_nd != 2) { + throw py::value_error("The array of pivot indices has ndim=" + + std::to_string(ipiv_array_nd) + + ", but a 2-dimensional array is expected."); + } + + // check compatibility of execution queue and allocation queue + if (!dpctl::utils::queues_are_compatible(exec_q, {a_array, ipiv_array})) { + throw py::value_error( + "Execution queue is not compatible with allocation queues"); + } + + auto const &overlap = dpctl::tensor::overlap::MemoryOverlap(); + if (overlap(a_array, ipiv_array)) { + throw py::value_error("The input array and the array of pivot indices " + "are overlapping segments of memory"); + } + + bool is_a_array_c_contig = a_array.is_c_contiguous(); + bool is_ipiv_array_c_contig = ipiv_array.is_c_contiguous(); + if (!is_a_array_c_contig) { + throw py::value_error("The input array " + "must be C-contiguous"); + } + if (!is_ipiv_array_c_contig) { + throw py::value_error("The array of pivot indices " + "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()); + + getrf_batch_impl_fn_ptr_t getrf_batch_fn = + getrf_batch_dispatch_vector[a_array_type_id]; + if (getrf_batch_fn == nullptr) { + throw py::value_error( + "No getrf_batch implementation defined for the provided type " + "of the input matrix."); + } + + auto ipiv_types = dpctl_td_ns::usm_ndarray_types(); + int ipiv_array_type_id = + ipiv_types.typenum_to_lookup_id(ipiv_array.get_typenum()); + + if (ipiv_array_type_id != static_cast(dpctl_td_ns::typenum_t::INT64)) { + throw py::value_error("The type of 'ipiv_array' must be int64."); + } + + char *a_array_data = a_array.get_data(); + const std::int64_t lda = std::max(1UL, n); + + char *ipiv_array_data = ipiv_array.get_data(); + std::int64_t *d_ipiv = reinterpret_cast(ipiv_array_data); + + std::vector host_task_events; + sycl::event getrf_batch_ev = getrf_batch_fn( + exec_q, n, a_array_data, lda, stride_a, d_ipiv, stride_ipiv, batch_size, + dev_info, host_task_events, depends); + + sycl::event args_ev = dpctl::utils::keep_args_alive( + exec_q, {a_array, ipiv_array}, host_task_events); + + return std::make_pair(args_ev, getrf_batch_ev); +} + +template +struct GetrfBatchContigFactory +{ + fnT get() + { + if constexpr (types::GetrfBatchTypePairSupportFactory::is_defined) { + return getrf_batch_impl; + } + else { + return nullptr; + } + } +}; + +void init_getrf_batch_dispatch_vector(void) +{ + dpctl_td_ns::DispatchVectorBuilder + contig; + contig.populate_dispatch_vector(getrf_batch_dispatch_vector); +} +} // namespace lapack +} // namespace ext +} // namespace backend +} // namespace dpnp diff --git a/dpnp/backend/extensions/lapack/lapack_py.cpp b/dpnp/backend/extensions/lapack/lapack_py.cpp index c0765be7509d..3c6ba80575da 100644 --- a/dpnp/backend/extensions/lapack/lapack_py.cpp +++ b/dpnp/backend/extensions/lapack/lapack_py.cpp @@ -31,6 +31,7 @@ #include #include "gesv.hpp" +#include "getrf.hpp" #include "heevd.hpp" #include "linalg_exceptions.hpp" #include "syevd.hpp" @@ -42,6 +43,8 @@ namespace py = pybind11; 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_syevd_dispatch_vector(); } @@ -68,6 +71,20 @@ PYBIND11_MODULE(_lapack_impl, m) py::arg("sycl_queue"), py::arg("coeff_matrix"), py::arg("dependent_vals"), py::arg("depends") = py::list()); + m.def("_getrf", &lapack_ext::getrf, + "Call `getrf` from OneMKL LAPACK library to return " + "the LU factorization of a general n x n matrix", + py::arg("sycl_queue"), py::arg("a_array"), py::arg("ipiv_array"), + py::arg("dev_info"), py::arg("depends") = py::list()); + + m.def("_getrf_batch", &lapack_ext::getrf_batch, + "Call `getrf_batch` from OneMKL LAPACK library to return " + "the LU factorization of a batch of general n x n matrices", + py::arg("sycl_queue"), py::arg("a_array"), py::arg("ipiv_array"), + py::arg("dev_info_array"), py::arg("n"), py::arg("stride_a"), + py::arg("stride_ipiv"), py::arg("batch_size"), + py::arg("depends") = py::list()); + m.def("_heevd", &lapack_ext::heevd, "Call `heevd` from OneMKL LAPACK library to return " "the eigenvalues and eigenvectors of a complex Hermitian matrix", diff --git a/dpnp/backend/extensions/lapack/types_matrix.hpp b/dpnp/backend/extensions/lapack/types_matrix.hpp index 60521cb75a3c..59198e975c8a 100644 --- a/dpnp/backend/extensions/lapack/types_matrix.hpp +++ b/dpnp/backend/extensions/lapack/types_matrix.hpp @@ -69,6 +69,59 @@ struct GesvTypePairSupportFactory // 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::getrf + * function. + * + * @tparam T Type of array containing input matrix, + * as well as the output array for storing the LU factorization. + */ +template +struct GetrfTypePairSupportFactory +{ + 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::getrf_batch + * function. + * + * @tparam T Type of array containing input matrix, + * as well as the output array for storing the LU factorization. + */ +template +struct GetrfBatchTypePairSupportFactory +{ + 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::heevd diff --git a/dpnp/backend/include/dpnp_iface_fptr.hpp b/dpnp/backend/include/dpnp_iface_fptr.hpp index aaef5c78627b..58e3fc106340 100644 --- a/dpnp/backend/include/dpnp_iface_fptr.hpp +++ b/dpnp/backend/include/dpnp_iface_fptr.hpp @@ -121,8 +121,6 @@ enum class DPNPFuncName : size_t DPNP_FN_DEGREES_EXT, /**< Used in numpy.degrees() impl, requires extra parameters */ DPNP_FN_DET, /**< Used in numpy.linalg.det() impl */ - DPNP_FN_DET_EXT, /**< Used in numpy.linalg.det() impl, requires extra - parameters */ DPNP_FN_DIAG, /**< Used in numpy.diag() impl */ DPNP_FN_DIAG_INDICES, /**< Used in numpy.diag_indices() impl */ DPNP_FN_DIAG_INDICES_EXT, /**< Used in numpy.diag_indices() impl, requires diff --git a/dpnp/backend/kernels/dpnp_krnl_indexing.cpp b/dpnp/backend/kernels/dpnp_krnl_indexing.cpp index 7dc35fb5a803..0889358d989f 100644 --- a/dpnp/backend/kernels/dpnp_krnl_indexing.cpp +++ b/dpnp/backend/kernels/dpnp_krnl_indexing.cpp @@ -935,6 +935,10 @@ void func_map_init_indexing_func(func_map_t &fmap) eft_FLT, (void *)dpnp_diagonal_ext_c}; fmap[DPNPFuncName::DPNP_FN_DIAGONAL_EXT][eft_DBL][eft_DBL] = { eft_DBL, (void *)dpnp_diagonal_ext_c}; + fmap[DPNPFuncName::DPNP_FN_DIAGONAL_EXT][eft_C64][eft_C64] = { + eft_C64, (void *)dpnp_diagonal_ext_c>}; + fmap[DPNPFuncName::DPNP_FN_DIAGONAL_EXT][eft_C128][eft_C128] = { + eft_C128, (void *)dpnp_diagonal_ext_c>}; fmap[DPNPFuncName::DPNP_FN_FILL_DIAGONAL][eft_INT][eft_INT] = { eft_INT, (void *)dpnp_fill_diagonal_default_c}; diff --git a/dpnp/backend/kernels/dpnp_krnl_linalg.cpp b/dpnp/backend/kernels/dpnp_krnl_linalg.cpp index 5e78a7cda176..8e970884a184 100644 --- a/dpnp/backend/kernels/dpnp_krnl_linalg.cpp +++ b/dpnp/backend/kernels/dpnp_krnl_linalg.cpp @@ -874,15 +874,6 @@ void func_map_init_linalg_func(func_map_t &fmap) fmap[DPNPFuncName::DPNP_FN_DET][eft_DBL][eft_DBL] = { eft_DBL, (void *)dpnp_det_default_c}; - fmap[DPNPFuncName::DPNP_FN_DET_EXT][eft_INT][eft_INT] = { - eft_INT, (void *)dpnp_det_ext_c}; - fmap[DPNPFuncName::DPNP_FN_DET_EXT][eft_LNG][eft_LNG] = { - eft_LNG, (void *)dpnp_det_ext_c}; - fmap[DPNPFuncName::DPNP_FN_DET_EXT][eft_FLT][eft_FLT] = { - eft_FLT, (void *)dpnp_det_ext_c}; - fmap[DPNPFuncName::DPNP_FN_DET_EXT][eft_DBL][eft_DBL] = { - eft_DBL, (void *)dpnp_det_ext_c}; - fmap[DPNPFuncName::DPNP_FN_INV][eft_INT][eft_INT] = { eft_DBL, (void *)dpnp_inv_default_c}; fmap[DPNPFuncName::DPNP_FN_INV][eft_LNG][eft_LNG] = { diff --git a/dpnp/dpnp_algo/dpnp_algo.pxd b/dpnp/dpnp_algo/dpnp_algo.pxd index 69aa29bf7172..cd3694bc0040 100644 --- a/dpnp/dpnp_algo/dpnp_algo.pxd +++ b/dpnp/dpnp_algo/dpnp_algo.pxd @@ -54,8 +54,6 @@ cdef extern from "dpnp_iface_fptr.hpp" namespace "DPNPFuncName": # need this na DPNP_FN_CUMSUM_EXT DPNP_FN_DEGREES DPNP_FN_DEGREES_EXT - DPNP_FN_DET - DPNP_FN_DET_EXT DPNP_FN_DIAG_INDICES DPNP_FN_DIAG_INDICES_EXT DPNP_FN_DIAGONAL diff --git a/dpnp/linalg/dpnp_algo_linalg.pyx b/dpnp/linalg/dpnp_algo_linalg.pyx index c86b869acd3c..ff5ea6253db6 100644 --- a/dpnp/linalg/dpnp_algo_linalg.pyx +++ b/dpnp/linalg/dpnp_algo_linalg.pyx @@ -47,7 +47,6 @@ cimport dpnp.dpnp_utils as utils __all__ = [ "dpnp_cholesky", "dpnp_cond", - "dpnp_det", "dpnp_eig", "dpnp_eigvals", "dpnp_inv", @@ -140,47 +139,6 @@ cpdef object dpnp_cond(object input, object p): return ret -cpdef utils.dpnp_descriptor dpnp_det(utils.dpnp_descriptor input): - cdef shape_type_c input_shape = input.shape - cdef size_t n = input.shape[-1] - cdef shape_type_c result_shape = (1,) - if input.ndim != 2: - result_shape = tuple((list(input.shape))[:-2]) - - cdef DPNPFuncType param1_type = dpnp_dtype_to_DPNPFuncType(input.dtype) - - cdef DPNPFuncData kernel_data = get_dpnp_function_ptr(DPNP_FN_DET_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(result_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_func_ptr_t func = kernel_data.ptr - - cdef c_dpctl.DPCTLSyclEventRef event_ref = func(q_ref, - input.get_data(), - result.get_data(), - input_shape.data(), - input.ndim, - NULL) # dep_events_ref - - with nogil: c_dpctl.DPCTLEvent_WaitAndThrow(event_ref) - c_dpctl.DPCTLEvent_Delete(event_ref) - - return result - - cpdef tuple dpnp_eig(utils.dpnp_descriptor x1): cdef shape_type_c x1_shape = x1.shape diff --git a/dpnp/linalg/dpnp_iface_linalg.py b/dpnp/linalg/dpnp_iface_linalg.py index c9c9c855728a..b7a56c4c7db2 100644 --- a/dpnp/linalg/dpnp_iface_linalg.py +++ b/dpnp/linalg/dpnp_iface_linalg.py @@ -50,7 +50,9 @@ from .dpnp_utils_linalg import ( check_stacked_2d, check_stacked_square, + dpnp_det, dpnp_eigh, + dpnp_slogdet, dpnp_solve, ) @@ -69,6 +71,7 @@ "qr", "solve", "svd", + "slogdet", ] @@ -148,32 +151,50 @@ def cond(input, p=None): return call_origin(numpy.linalg.cond, input, p) -def det(input): +def det(a): """ Compute the determinant of an array. + For full documentation refer to :obj:`numpy.linalg.det`. + Parameters ---------- - input : (..., M, M) array_like + a : (..., M, M) {dpnp.ndarray, usm_ndarray} Input array to compute determinants for. Returns ------- - det : (...) array_like - Determinant of `input`. - """ + det : (...) dpnp.ndarray + Determinant of `a`. - x1_desc = dpnp.get_dpnp_descriptor(input, copy_when_nondefault_queue=False) - if x1_desc: - if x1_desc.ndim < 2: - pass - elif x1_desc.shape[-1] == x1_desc.shape[-2]: - result_obj = dpnp_det(x1_desc).get_pyobj() - result = dpnp.convert_single_elem_array_to_scalar(result_obj) + See Also + -------- + :obj:`dpnp.linalg.slogdet` : Returns sign and logarithm of the determinant of an array. - return result + Examples + -------- + The determinant of a 2-D array [[a, b], [c, d]] is ad - bc: - return call_origin(numpy.linalg.det, input) + >>> import dpnp as dp + >>> a = dp.array([[1, 2], [3, 4]]) + >>> dp.linalg.det(a) + array(-2.) + + Computing determinants for a stack of matrices: + + >>> a = dp.array([ [[1, 2], [3, 4]], [[1, 2], [2, 1]], [[1, 3], [3, 1]] ]) + >>> a.shape + (3, 2, 2) + >>> dp.linalg.det(a) + array([-2., -3., -8.]) + + """ + + dpnp.check_supported_arrays_type(a) + check_stacked_2d(a) + check_stacked_square(a) + + return dpnp_det(a) def eig(x1): @@ -633,3 +654,59 @@ def svd(x1, full_matrices=True, compute_uv=True, hermitian=False): return call_origin( numpy.linalg.svd, x1, full_matrices, compute_uv, hermitian ) + + +def slogdet(a): + """ + Compute the sign and (natural) logarithm of the determinant of an array. + + For full documentation refer to :obj:`numpy.linalg.slogdet`. + + Parameters + ---------- + a : (..., M, M) {dpnp.ndarray, usm_ndarray} + Input array, has to be a square 2-D array. + + Returns + ------- + sign : (...) dpnp.ndarray + A number representing the sign of the determinant. For a real matrix, + this is 1, 0, or -1. For a complex matrix, this is a complex number + with absolute value 1 (i.e., it is on the unit circle), or else 0. + logabsdet : (...) dpnp.ndarray + The natural log of the absolute value of the determinant. + + See Also + -------- + :obj:`dpnp.det` : Returns the determinant of an array. + + Examples + -------- + The determinant of a 2-D array [[a, b], [c, d]] is ad - bc: + + >>> import dpnp as dp + >>> a = dp.array([[1, 2], [3, 4]]) + >>> (sign, logabsdet) = dp.linalg.slogdet(a) + >>> (sign, logabsdet) + (array(-1.), array(0.69314718)) + >>> sign * dp.exp(logabsdet) + array(-2.) + + Computing log-determinants for a stack of matrices: + + >>> a = dp.array([ [[1, 2], [3, 4]], [[1, 2], [2, 1]], [[1, 3], [3, 1]] ]) + >>> a.shape + (3, 2, 2) + >>> sign, logabsdet = dp.linalg.slogdet(a) + >>> (sign, logabsdet) + (array([-1., -1., -1.]), array([0.69314718, 1.09861229, 2.07944154])) + >>> sign * dp.exp(logabsdet) + array([-2., -3., -8.]) + + """ + + dpnp.check_supported_arrays_type(a) + check_stacked_2d(a) + check_stacked_square(a) + + return dpnp_slogdet(a) diff --git a/dpnp/linalg/dpnp_utils_linalg.py b/dpnp/linalg/dpnp_utils_linalg.py index a3a2802072c9..bd44b9742534 100644 --- a/dpnp/linalg/dpnp_utils_linalg.py +++ b/dpnp/linalg/dpnp_utils_linalg.py @@ -37,13 +37,89 @@ __all__ = [ "check_stacked_2d", "check_stacked_square", + "dpnp_det", "dpnp_eigh", + "dpnp_slogdet", "dpnp_solve", ] _jobz = {"N": 0, "V": 1} _upper_lower = {"U": 0, "L": 1} +_real_types_map = { + "float32": "float32", # single : single + "float64": "float64", # double : double + "complex64": "float32", # csingle : csingle + "complex128": "float64", # cdouble : cdouble +} + + +def _calculate_determinant_sign(ipiv, diag, res_type, n): + """ + Calculate the sign of the determinant based on row exchanges and diagonal values. + + Parameters + ----------- + ipiv : {dpnp.ndarray, usm_ndarray} + The pivot indices from LU decomposition. + diag : {dpnp.ndarray, usm_ndarray} + The diagonal elements of the LU decomposition matrix. + res_type : dpnp.dtype + The common data type for linalg operations. + n : int + The size of the last two dimensions of the array. + + Returns + ------- + sign : {dpnp_array, usm_ndarray} + The sign of the determinant. + + """ + + # Checks for row exchanges in LU decomposition affecting determinant sign. + ipiv_diff = ipiv != dpnp.arange( + 1, n + 1, usm_type=ipiv.usm_type, sycl_queue=ipiv.sycl_queue + ) + + # Counts row exchanges from 'ipiv_diff'. + non_zero = dpnp.count_nonzero(ipiv_diff, axis=-1) + + # For floating types, adds count of negative diagonal elements + # to determine determinant sign. + if dpnp.issubdtype(res_type, dpnp.floating): + non_zero += dpnp.count_nonzero(diag < 0, axis=-1) + + sign = (non_zero % 2) * -2 + 1 + + # For complex types, compute sign from the phase of diagonal elements. + if dpnp.issubdtype(res_type, dpnp.complexfloating): + sign = sign * dpnp.prod(diag / dpnp.abs(diag), axis=-1) + + return sign.astype(res_type) + + +def _real_type(dtype, device=None): + """ + Returns the real data type corresponding to a given dpnp data type. + + Parameters + ---------- + dtype : dpnp.dtype + The dtype for which to find the corresponding real data type. + device : {None, string, SyclDevice, SyclQueue}, optional + An array API concept of device where an array of default floating type might be created. + + Returns + ------- + out : str + The name of the real data type. + + """ + + default = dpnp.default_float_type(device) + real_type = _real_types_map.get(dtype.name, default) + return dpnp.dtype(real_type) + def check_stacked_2d(*arrays): """ @@ -53,7 +129,7 @@ def check_stacked_2d(*arrays): Parameters ---------- - arrays : {dpnp_array, usm_ndarray} + arrays : {dpnp.ndarray, usm_ndarray} A sequence of input arrays to check for dimensionality. Returns @@ -92,7 +168,7 @@ def check_stacked_square(*arrays): Parameters ---------- - arrays : {dpnp_array, usm_ndarray} + arrays : {dpnp.ndarray, usm_ndarray} A sequence of input arrays to check for square matrix shape. Returns @@ -117,8 +193,6 @@ def check_stacked_square(*arrays): def _common_type(*arrays): """ - _common_type(*arrays) - Common type for linear algebra operations. This function determines the common data type for linalg operations. @@ -129,12 +203,15 @@ def _common_type(*arrays): - The default floating-point data type is determined by the capabilities of the device on which `arrays` are created, as indicated by `dpnp.default_float_type()`. - Args: - *arrays (dpnp.ndarray): Input arrays. - - Returns: - dtype_common (dtype): The common data type for linalg operations. + Parameters + ---------- + arrays : {dpnp.ndarray, usm_ndarray} + A sequence of input arrays. + Returns + ------- + dtype_common : dpnp.dtype + The common data type for linalg operations. This returned value is applicable both as the precision to be used in linalg calls and as the dtype of (possibly complex) output(s). @@ -150,30 +227,269 @@ def _common_type(*arrays): def _common_inexact_type(default_dtype, *dtypes): """ - _common_inexact_type(default_dtype, *dtypes) - Determines the common 'inexact' data type for linear algebra operations. This function selects an 'inexact' data type appropriate for the device's capabilities. It defaults to `default_dtype` when provided types are not 'inexact'. - Args: - default_dtype: The default data type. This is determined by the capabilities of + Parameters + ---------- + default_dtype : dpnp.dtype + The default data type. This is determined by the capabilities of the device and is used when none of the provided types are 'inexact'. *dtypes: A variable number of data types to be evaluated to find the common 'inexact' type. - Returns: - dpnp.result_type (dtype) : The resultant 'inexact' data type for linalg operations, + Returns + ------- + dpnp.result_type : dpnp.dtype + The resultant 'inexact' data type for linalg operations, ensuring computational compatibility. """ + inexact_dtypes = [ dt if issubdtype(dt, dpnp.inexact) else default_dtype for dt in dtypes ] return dpnp.result_type(*inexact_dtypes) +def _lu_factor(a, res_type): + """ + Compute pivoted LU decomposition. + + Decompose a given batch of square matrices. Inputs and outputs are + transposed. + + Parameters + ---------- + a : (..., M, M) {dpnp.ndarray, usm_ndarray} + Input array containing the matrices to be decomposed. + res_type : dpnp.dtype + Specifies the data type of the result. + Acceptable data types are float32, float64, complex64, or complex128. + + Returns + ------- + tuple: + lu_t : (..., N, N) {dpnp.ndarray, usm_ndarray} + Combined 'L' and 'U' matrices from LU decomposition + excluding the diagonal of 'L'. + piv : (..., N) {dpnp.ndarray, usm_ndarray} + 1-origin pivot indices indicating row permutations during decomposition. + dev_info : (...) {dpnp.ndarray, usm_ndarray} + Information on `getrf` or `getrf_batch` computation success (0 for success). + + """ + + n = a.shape[-2] + + a_sycl_queue = a.sycl_queue + a_usm_type = a.usm_type + + # TODO: Find out at which array sizes the best performance is obtained + # getrf_batch implementation shows slow results with large arrays on GPU. + # Use getrf_batch only on CPU. + # On GPU call getrf for each two-dimensional array by loop + use_batch = a.sycl_device.has_aspect_cpu + + if a.ndim > 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) + + if use_batch: + # `a` must be copied because getrf_batch destroys the input matrix + a_h = dpnp.empty_like(a, order="C", dtype=res_type) + ipiv_h = dpnp.empty( + (batch_size, n), + dtype=dpnp.int64, + order="C", + usm_type=a_usm_type, + sycl_queue=a_sycl_queue, + ) + dev_info_h = [0] * batch_size + + 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 + ) + + ipiv_stride = n + a_stride = a_h.strides[0] + + # Call the LAPACK extension function _getrf_batch + # to perform LU decomposition of a batch of general matrices + ht_lapack_ev, _ = li._getrf_batch( + a_sycl_queue, + a_h.get_array(), + ipiv_h.get_array(), + dev_info_h, + n, + a_stride, + ipiv_stride, + batch_size, + [a_copy_ev], + ) + + ht_lapack_ev.wait() + a_ht_copy_ev.wait() + + dev_info_array = dpnp.array( + dev_info_h, usm_type=a_usm_type, sycl_queue=a_sycl_queue + ) + + # Reshape the results back to their original shape + a_h = a_h.reshape(orig_shape) + ipiv_h = ipiv_h.reshape(orig_shape[:-1]) + dev_info_array = dev_info_array.reshape(orig_shape[:-2]) + + return (a_h, ipiv_h, dev_info_array) + + else: + # Initialize lists for storing arrays and events for each batch + a_vecs = [None] * batch_size + ipiv_vecs = [None] * batch_size + dev_info_vecs = [None] * batch_size + a_ht_copy_ev = [None] * batch_size + ht_lapack_ev = [None] * batch_size + + # Process each batch + for i in range(batch_size): + # Copy each 2D slice to a new array as getrf destroys the input matrix + a_vecs[i] = dpnp.empty_like(a[i], order="C", dtype=res_type) + ( + a_ht_copy_ev[i], + a_copy_ev, + ) = ti._copy_usm_ndarray_into_usm_ndarray( + src=a_usm_arr[i], + dst=a_vecs[i].get_array(), + sycl_queue=a_sycl_queue, + ) + ipiv_vecs[i] = dpnp.empty( + (n,), + dtype=dpnp.int64, + order="C", + usm_type=a_usm_type, + sycl_queue=a_sycl_queue, + ) + dev_info_vecs[i] = [0] + + # Call the LAPACK extension function _getrf + # to perform LU decomposition on each batch in 'a_vecs[i]' + ht_lapack_ev[i], _ = li._getrf( + a_sycl_queue, + a_vecs[i].get_array(), + ipiv_vecs[i].get_array(), + dev_info_vecs[i], + [a_copy_ev], + ) + + for i in range(batch_size): + ht_lapack_ev[i].wait() + a_ht_copy_ev[i].wait() + + # Reshape the results back to their original shape + out_a = dpnp.array(a_vecs, order="C").reshape(orig_shape) + out_ipiv = dpnp.array(ipiv_vecs).reshape(orig_shape[:-1]) + out_dev_info = dpnp.array( + dev_info_vecs, usm_type=a_usm_type, sycl_queue=a_sycl_queue + ).reshape(orig_shape[:-2]) + + return (out_a, out_ipiv, out_dev_info) + + else: + a_usm_arr = dpnp.get_usm_ndarray(a) + + # `a` must be copied because getrf destroys the input matrix + a_h = dpnp.empty_like(a, order="C", dtype=res_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 + ) + + ipiv_h = dpnp.empty( + n, + dtype=dpnp.int64, + order="C", + usm_type=a_usm_type, + sycl_queue=a_sycl_queue, + ) + dev_info_h = [0] + + # Call the LAPACK extension function _getrf + # to perform LU decomposition on the input matrix + ht_lapack_ev, _ = li._getrf( + a_sycl_queue, + a_h.get_array(), + ipiv_h.get_array(), + dev_info_h, + [a_copy_ev], + ) + + ht_lapack_ev.wait() + a_ht_copy_ev.wait() + + dev_info_array = dpnp.array( + dev_info_h, usm_type=a_usm_type, sycl_queue=a_sycl_queue + ) + + # Return a tuple containing the factorized matrix 'a_h', + # pivot indices 'ipiv_h' + # and the status 'dev_info_h' from the LAPACK getrf call + return (a_h, ipiv_h, dev_info_array) + + +def dpnp_det(a): + """ + dpnp_det(a) + + Returns the determinant of `a` array. + + """ + + a_usm_type = a.usm_type + a_sycl_queue = a.sycl_queue + + res_type = _common_type(a) + + a_shape = a.shape + shape = a_shape[:-2] + n = a_shape[-2] + + if a.size == 0: + # empty batch (result is empty, too) or empty matrices det([[]]) == 1 + det = dpnp.ones( + shape, + dtype=res_type, + usm_type=a_usm_type, + sycl_queue=a_sycl_queue, + ) + return det + + lu, ipiv, dev_info = _lu_factor(a, res_type) + + # Transposing 'lu' to swap the last two axes for compatibility + # with 'dpnp.diagonal' as it does not support 'axis1' and 'axis2' arguments. + # TODO: Replace with 'dpnp.diagonal(lu, axis1=-2, axis2=-1)' when supported. + lu_transposed = lu.transpose(-2, -1, *range(lu.ndim - 2)) + diag = dpnp.diagonal(lu_transposed) + + det = dpnp.prod(dpnp.abs(diag), axis=-1) + + sign = _calculate_determinant_sign(ipiv, diag, res_type, n) + + det = sign * det + det = det.astype(res_type, copy=False) + singular = dev_info > 0 + det = dpnp.where(singular, res_type.type(0), det) + + return det.reshape(shape) + + def dpnp_eigh(a, UPLO): """ dpnp_eigh(a, UPLO) @@ -432,3 +748,54 @@ def dpnp_solve(a, b): a_ht_copy_ev.wait() return b_f + + +def dpnp_slogdet(a): + """ + dpnp_slogdet(a) + + Returns sign and logarithm of the determinant of `a` array. + + """ + + a_usm_type = a.usm_type + a_sycl_queue = a.sycl_queue + + res_type = _common_type(a) + logdet_dtype = _real_type(res_type) + + a_shape = a.shape + shape = a_shape[:-2] + n = a_shape[-2] + + if a.size == 0: + # empty batch (result is empty, too) or empty matrices det([[]]) == 1 + sign = dpnp.ones( + shape, dtype=res_type, usm_type=a_usm_type, sycl_queue=a_sycl_queue + ) + logdet = dpnp.zeros( + shape, + dtype=logdet_dtype, + usm_type=a_usm_type, + sycl_queue=a_sycl_queue, + ) + return sign, logdet + + lu, ipiv, dev_info = _lu_factor(a, res_type) + + # Transposing 'lu' to swap the last two axes for compatibility + # with 'dpnp.diagonal' as it does not support 'axis1' and 'axis2' arguments. + # TODO: Replace with 'dpnp.diagonal(lu, axis1=-2, axis2=-1)' when supported. + lu_transposed = lu.transpose(-2, -1, *range(lu.ndim - 2)) + diag = dpnp.diagonal(lu_transposed) + + logdet = dpnp.log(dpnp.abs(diag)).sum(axis=-1) + + sign = _calculate_determinant_sign(ipiv, diag, res_type, n) + + logdet = logdet.astype(logdet_dtype, copy=False) + singular = dev_info > 0 + return ( + dpnp.where(singular, res_type.type(0), sign).reshape(shape), + dpnp.where(singular, logdet_dtype.type("-inf"), logdet).reshape(shape), + ) diff --git a/tests/test_indexing.py b/tests/test_indexing.py index 4d8229e53ce7..3cb9d4298d61 100644 --- a/tests/test_indexing.py +++ b/tests/test_indexing.py @@ -261,6 +261,7 @@ def test_choose(): assert_array_equal(expected, result) +@pytest.mark.parametrize("arr_dtype", get_all_dtypes(no_bool=True)) @pytest.mark.parametrize("offset", [0, 1], ids=["0", "1"]) @pytest.mark.parametrize( "array", @@ -296,8 +297,8 @@ def test_choose(): "[[[[1, 2, 3], [4, 5, 6]], [[7, 8, 9], [10, 11, 12]]], [[[13, 14, 15], [16, 17, 18]], [[19, 20, 21], [22, 23, 24]]]]", ], ) -def test_diagonal(array, offset): - a = numpy.array(array) +def test_diagonal(array, arr_dtype, offset): + a = numpy.array(array, dtype=arr_dtype) ia = dpnp.array(a) expected = numpy.diagonal(a, offset) result = dpnp.diagonal(ia, offset) diff --git a/tests/test_linalg.py b/tests/test_linalg.py index 6c8c2a5f04e3..cb780abd9253 100644 --- a/tests/test_linalg.py +++ b/tests/test_linalg.py @@ -105,46 +105,118 @@ def test_cond(arr, p): assert_array_equal(expected, result) -@pytest.mark.parametrize( - "array", - [ - [[0, 0], [0, 0]], - [[1, 2], [1, 2]], - [[1, 2], [3, 4]], - [[[1, 2], [3, 4]], [[1, 2], [2, 1]], [[1, 3], [3, 1]]], +class TestDet: + # TODO: Remove the use of fixture for test_det + # when dpnp.prod() will support complex dtypes on Gen9 + @pytest.mark.usefixtures("allow_fall_back_on_numpy") + @pytest.mark.parametrize( + "array", [ - [[[1, 2], [3, 4]], [[1, 2], [2, 1]]], - [[[1, 3], [3, 1]], [[0, 1], [1, 3]]], + [[1, 2], [3, 4]], + [[[1, 2], [3, 4]], [[1, 2], [2, 1]], [[1, 3], [3, 1]]], + [ + [[[1, 2], [3, 4]], [[1, 2], [2, 1]]], + [[[1, 3], [3, 1]], [[0, 1], [1, 3]]], + ], ], - ], - ids=[ - "[[0, 0], [0, 0]]", - "[[1, 2], [1, 2]]", - "[[1, 2], [3, 4]]", - "[[[1, 2], [3, 4]], [[1, 2], [2, 1]], [[1, 3], [3, 1]]]", - "[[[[1, 2], [3, 4]], [[1, 2], [2, 1]]], [[[1, 3], [3, 1]], [[0, 1], [1, 3]]]]", - ], -) -def test_det(array): - a = numpy.array(array) - ia = inp.array(a) - result = inp.linalg.det(ia) - expected = numpy.linalg.det(a) - assert_allclose(expected, result) + ids=[ + "2D_array", + "3D_array", + "4D_array", + ], + ) + @pytest.mark.parametrize("dtype", get_all_dtypes(no_bool=True)) + def test_det(self, array, dtype): + a = numpy.array(array, dtype=dtype) + ia = inp.array(a) + result = inp.linalg.det(ia) + expected = numpy.linalg.det(a) + assert_allclose(expected, result) + def test_det_strides(self): + a_np = numpy.array( + [ + [2, 3, 1, 4, 5], + [5, 6, 7, 8, 9], + [9, 7, 7, 2, 3], + [1, 4, 5, 1, 8], + [8, 9, 8, 5, 3], + ] + ) -@pytest.mark.usefixtures("allow_fall_back_on_numpy") -def test_det_empty(): - a = numpy.empty((0, 0, 2, 2), dtype=numpy.float32) - ia = inp.array(a) + a_dp = inp.array(a_np) + + # positive strides + expected = numpy.linalg.det(a_np[::2, ::2]) + result = inp.linalg.det(a_dp[::2, ::2]) + assert_allclose(expected, result, rtol=1e-3, atol=1e-4) + + # negative strides + expected = numpy.linalg.det(a_np[::-2, ::-2]) + result = inp.linalg.det(a_dp[::-2, ::-2]) + assert_allclose(expected, result, rtol=1e-3, atol=1e-4) + + def test_det_empty(self): + a = numpy.empty((0, 0, 2, 2), dtype=numpy.float32) + ia = inp.array(a) + + np_det = numpy.linalg.det(a) + dpnp_det = inp.linalg.det(ia) + + assert dpnp_det.dtype == np_det.dtype + assert dpnp_det.shape == np_det.shape + + assert_allclose(np_det, dpnp_det) + + @pytest.mark.parametrize( + "matrix", + [ + [[1, 2], [2, 4]], + [[0, 0], [0, 0]], + [[1, 1], [1, 1]], + [[2, 4], [1, 2]], + [[1, 2], [0, 0]], + [[1, 0], [2, 0]], + ], + ids=[ + "Linearly dependent rows", + "Zero matrix", + "Identical rows", + "Linearly dependent columns", + "Zero row", + "Zero column", + ], + ) + def test_det_singular_matrix(self, matrix): + a_np = numpy.array(matrix, dtype="float32") + a_dp = inp.array(a_np) - np_det = numpy.linalg.det(a) - dpnp_det = inp.linalg.det(ia) + expected = numpy.linalg.slogdet(a_np) + result = inp.linalg.slogdet(a_dp) - assert dpnp_det.dtype == np_det.dtype - assert dpnp_det.shape == np_det.shape + assert_allclose(expected, result, rtol=1e-3, atol=1e-4) + + # TODO: remove skipif when MKLD-16626 is resolved + # _getrf_batch does not raise an error with singular matrices. + # Skip running on cpu because dpnp uses _getrf_batch only on cpu. + @pytest.mark.skipif(is_cpu_device(), reason="MKLD-16626") + def test_det_singular_matrix_3D(self): + a_np = numpy.array( + [[[1, 2], [3, 4]], [[1, 2], [1, 2]], [[1, 3], [3, 1]]] + ) + a_dp = inp.array(a_np) - assert_allclose(np_det, dpnp_det) + expected = numpy.linalg.det(a_np) + result = inp.linalg.det(a_dp) + + assert_allclose(expected, result, rtol=1e-3, atol=1e-4) + + def test_det_errors(self): + a_dp = inp.array([[1, 2], [3, 5]], dtype="float32") + + # unsupported type + a_np = inp.asnumpy(a_dp) + assert_raises(TypeError, inp.linalg.det, a_np) @pytest.mark.parametrize("type", get_all_dtypes(no_bool=True, no_complex=True)) @@ -618,3 +690,115 @@ def test_solve_errors(self): assert_raises( inp.linalg.LinAlgError, inp.linalg.solve, a_dp_ndim_1, b_dp ) + + +class TestSlogdet: + # TODO: Remove the use of fixture for test_slogdet_2d and test_slogdet_3d + # when dpnp.prod() will support complex dtypes on Gen9 + @pytest.mark.usefixtures("allow_fall_back_on_numpy") + @pytest.mark.parametrize("dtype", get_all_dtypes(no_bool=True)) + def test_slogdet_2d(self, dtype): + a_np = numpy.array([[1, 2], [3, 4]], dtype=dtype) + a_dp = inp.array(a_np) + + sign_expected, logdet_expected = numpy.linalg.slogdet(a_np) + sign_result, logdet_result = inp.linalg.slogdet(a_dp) + + assert_allclose(sign_expected, sign_result) + assert_allclose(logdet_expected, logdet_result, rtol=1e-3, atol=1e-4) + + @pytest.mark.usefixtures("allow_fall_back_on_numpy") + @pytest.mark.parametrize("dtype", get_all_dtypes(no_bool=True)) + def test_slogdet_3d(self, dtype): + a_np = numpy.array( + [ + [[1, 2], [3, 4]], + [[1, 2], [2, 1]], + [[1, 3], [3, 1]], + ], + dtype=dtype, + ) + a_dp = inp.array(a_np) + + sign_expected, logdet_expected = numpy.linalg.slogdet(a_np) + sign_result, logdet_result = inp.linalg.slogdet(a_dp) + + assert_allclose(sign_expected, sign_result) + assert_allclose(logdet_expected, logdet_result, rtol=1e-3, atol=1e-4) + + def test_slogdet_strides(self): + a_np = numpy.array( + [ + [2, 3, 1, 4, 5], + [5, 6, 7, 8, 9], + [9, 7, 7, 2, 3], + [1, 4, 5, 1, 8], + [8, 9, 8, 5, 3], + ] + ) + + a_dp = inp.array(a_np) + + # positive strides + sign_expected, logdet_expected = numpy.linalg.slogdet(a_np[::2, ::2]) + sign_result, logdet_result = inp.linalg.slogdet(a_dp[::2, ::2]) + assert_allclose(sign_expected, sign_result) + assert_allclose(logdet_expected, logdet_result, rtol=1e-3, atol=1e-4) + + # negative strides + sign_expected, logdet_expected = numpy.linalg.slogdet(a_np[::-2, ::-2]) + sign_result, logdet_result = inp.linalg.slogdet(a_dp[::-2, ::-2]) + assert_allclose(sign_expected, sign_result) + assert_allclose(logdet_expected, logdet_result, rtol=1e-3, atol=1e-4) + + @pytest.mark.parametrize( + "matrix", + [ + [[1, 2], [2, 4]], + [[0, 0], [0, 0]], + [[1, 1], [1, 1]], + [[2, 4], [1, 2]], + [[1, 2], [0, 0]], + [[1, 0], [2, 0]], + ], + ids=[ + "Linearly dependent rows", + "Zero matrix", + "Identical rows", + "Linearly dependent columns", + "Zero row", + "Zero column", + ], + ) + def test_slogdet_singular_matrix(self, matrix): + a_np = numpy.array(matrix, dtype="float32") + a_dp = inp.array(a_np) + + sign_expected, logdet_expected = numpy.linalg.slogdet(a_np) + sign_result, logdet_result = inp.linalg.slogdet(a_dp) + + assert_allclose(sign_expected, sign_result) + assert_allclose(logdet_expected, logdet_result, rtol=1e-3, atol=1e-4) + + # TODO: remove skipif when MKLD-16626 is resolved + # _getrf_batch does not raise an error with singular matrices. + # Skip running on cpu because dpnp uses _getrf_batch only on cpu. + @pytest.mark.skipif(is_cpu_device(), reason="MKLD-16626") + def test_slogdet_singular_matrix_3D(self): + a_np = numpy.array( + [[[1, 2], [3, 4]], [[1, 2], [1, 2]], [[1, 3], [3, 1]]] + ) + a_dp = inp.array(a_np) + + sign_expected, logdet_expected = numpy.linalg.slogdet(a_np) + sign_result, logdet_result = inp.linalg.slogdet(a_dp) + + assert_allclose(sign_expected, sign_result) + assert_allclose(logdet_expected, logdet_result, rtol=1e-3, atol=1e-4) + + def test_slogdet_errors(self): + a_dp = inp.array([[1, 2], [3, 5]], dtype="float32") + + # unsupported type + a_np = inp.asnumpy(a_dp) + assert_raises(TypeError, inp.linalg.slogdet, a_np) diff --git a/tests/test_sycl_queue.py b/tests/test_sycl_queue.py index 53bc37f4eebd..de762ecef93d 100644 --- a/tests/test_sycl_queue.py +++ b/tests/test_sycl_queue.py @@ -1399,6 +1399,17 @@ def test_diff_scalar_append(device, kwargs): assert_sycl_queue_equal(result_queue, expected_queue) +@pytest.mark.parametrize( + "device", + valid_devices, + ids=[device.filter_string for device in valid_devices], +) +def test_clip(device): + x = dpnp.arange(10, device=device) + y = dpnp.clip(x, 3, 7) + assert_sycl_queue_equal(x.sycl_queue, y.sycl_queue) + + @pytest.mark.parametrize("func", ["take", "take_along_axis"]) @pytest.mark.parametrize( "device", @@ -1445,12 +1456,44 @@ def test_solve(device): assert_sycl_queue_equal(result_queue, dpnp_y.sycl_queue) +@pytest.mark.parametrize( + "shape, is_empty", + [ + ((2, 2), False), + ((3, 2, 2), False), + ((0, 0), True), + ((0, 2, 2), True), + ], + ids=[ + "(2, 2)", + "(3, 2, 2)", + "(0, 0)", + "(0, 2, 2)", + ], +) @pytest.mark.parametrize( "device", valid_devices, ids=[device.filter_string for device in valid_devices], ) -def test_clip(device): - x = dpnp.arange(10, device=device) - y = dpnp.clip(x, 3, 7) - assert_sycl_queue_equal(x.sycl_queue, y.sycl_queue) +def test_slogdet(shape, is_empty, device): + if is_empty: + numpy_x = numpy.empty(shape, dtype=dpnp.default_float_type(device)) + else: + count_elem = numpy.prod(shape) + numpy_x = numpy.arange( + 1, count_elem + 1, dtype=dpnp.default_float_type(device) + ).reshape(shape) + + dpnp_x = dpnp.array(numpy_x, device=device) + + sign_result, logdet_result = dpnp.linalg.slogdet(dpnp_x) + sign_expected, logdet_expected = numpy.linalg.slogdet(numpy_x) + assert_allclose(logdet_expected, logdet_result, rtol=1e-3, atol=1e-4) + assert_allclose(sign_expected, sign_result) + + sign_queue = sign_result.sycl_queue + logdet_queue = logdet_result.sycl_queue + + assert_sycl_queue_equal(sign_queue, dpnp_x.sycl_queue) + assert_sycl_queue_equal(logdet_queue, dpnp_x.sycl_queue) diff --git a/tests/test_usm_type.py b/tests/test_usm_type.py index 144b58c9a6ba..6b46dcceb1e8 100644 --- a/tests/test_usm_type.py +++ b/tests/test_usm_type.py @@ -525,6 +525,13 @@ 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("usm_type", list_of_usm_types, ids=list_of_usm_types) +def test_clip(usm_type): + x = dp.arange(10, usm_type=usm_type) + y = dp.clip(x, 2, 7) + assert x.usm_type == y.usm_type + + @pytest.mark.parametrize( "usm_type_matrix", list_of_usm_types, ids=list_of_usm_types ) @@ -564,7 +571,61 @@ def test_solve(matrix, vector, usm_type_matrix, usm_type_vector): @pytest.mark.parametrize("usm_type", list_of_usm_types, ids=list_of_usm_types) -def test_clip(usm_type): - x = dp.arange(10, usm_type=usm_type) - y = dp.clip(x, 2, 7) - assert x.usm_type == y.usm_type +@pytest.mark.parametrize( + "shape, is_empty", + [ + ((2, 2), False), + ((3, 2, 2), False), + ((0, 0), True), + ((0, 2, 2), True), + ], + ids=[ + "(2, 2)", + "(3, 2, 2)", + "(0, 0)", + "(0, 2, 2)", + ], +) +def test_slogdet(shape, is_empty, usm_type): + if is_empty: + x = dp.empty(shape, dtype=dp.default_float_type(), usm_type=usm_type) + else: + count_elem = numpy.prod(shape) + x = dp.arange( + 1, count_elem + 1, dtype=dp.default_float_type(), usm_type=usm_type + ).reshape(shape) + + sign, logdet = dp.linalg.slogdet(x) + + assert x.usm_type == sign.usm_type + assert x.usm_type == logdet.usm_type + + +@pytest.mark.parametrize("usm_type", list_of_usm_types, ids=list_of_usm_types) +@pytest.mark.parametrize( + "shape, is_empty", + [ + ((2, 2), False), + ((3, 2, 2), False), + ((0, 0), True), + ((0, 2, 2), True), + ], + ids=[ + "(2, 2)", + "(3, 2, 2)", + "(0, 0)", + "(0, 2, 2)", + ], +) +def test_det(shape, is_empty, usm_type): + if is_empty: + x = dp.empty(shape, dtype=dp.default_float_type(), usm_type=usm_type) + else: + count_elem = numpy.prod(shape) + x = dp.arange( + 1, count_elem + 1, dtype=dp.default_float_type(), usm_type=usm_type + ).reshape(shape) + + det = dp.linalg.det(x) + + assert x.usm_type == det.usm_type diff --git a/tests/third_party/cupy/linalg_tests/test_norms.py b/tests/third_party/cupy/linalg_tests/test_norms.py new file mode 100644 index 000000000000..2ed49d16057a --- /dev/null +++ b/tests/third_party/cupy/linalg_tests/test_norms.py @@ -0,0 +1,136 @@ +import unittest + +import numpy +import pytest + +import dpnp as cupy +from tests.helper import is_cpu_device +from tests.third_party.cupy import testing + + +# TODO: Remove the use of fixture for all tests in this file +# when dpnp.prod() will support complex dtypes on Gen9 +@pytest.mark.usefixtures("allow_fall_back_on_numpy") +class TestDet(unittest.TestCase): + @testing.for_dtypes("fdFD") + @testing.numpy_cupy_allclose(rtol=1e-3, atol=1e-4) + def test_det(self, xp, dtype): + a = testing.shaped_arange((2, 2), xp, dtype) + 1 + return xp.linalg.det(a) + + @testing.for_dtypes("fdFD") + @testing.numpy_cupy_allclose(rtol=1e-3, atol=1e-4) + def test_det_3(self, xp, dtype): + a = testing.shaped_arange((2, 2, 2), xp, dtype) + 1 + return xp.linalg.det(a) + + @testing.for_dtypes("fdFD") + @testing.numpy_cupy_allclose(rtol=1e-3, atol=1e-4) + def test_det_4(self, xp, dtype): + a = testing.shaped_arange((2, 2, 2, 2), xp, dtype) + 1 + return xp.linalg.det(a) + + @testing.for_dtypes("fdFD") + @testing.numpy_cupy_allclose(rtol=1e-3, atol=1e-4) + def test_det_empty_batch(self, xp, dtype): + a = xp.empty((2, 0, 3, 3), dtype=dtype) + return xp.linalg.det(a) + + @testing.for_dtypes("fdFD") + @testing.numpy_cupy_allclose(rtol=1e-3, atol=1e-4) + def test_det_empty_matrix(self, xp, dtype): + a = xp.empty((0, 0), dtype=dtype) + return xp.linalg.det(a) + + @testing.for_dtypes("fdFD") + @testing.numpy_cupy_allclose(rtol=1e-3, atol=1e-4) + def test_det_empty_matrices(self, xp, dtype): + a = xp.empty((2, 3, 0, 0), dtype=dtype) + return xp.linalg.det(a) + + @testing.for_dtypes("fdFD") + def test_det_different_last_two_dims(self, dtype): + for xp in (numpy, cupy): + a = testing.shaped_arange((2, 3, 2), xp, dtype) + with pytest.raises(xp.linalg.LinAlgError): + xp.linalg.det(a) + + @testing.for_dtypes("fdFD") + def test_det_different_last_two_dims_empty_batch(self, dtype): + for xp in (numpy, cupy): + a = xp.empty((0, 3, 2), dtype=dtype) + with pytest.raises(xp.linalg.LinAlgError): + xp.linalg.det(a) + + @testing.for_dtypes("fdFD") + def test_det_one_dim(self, dtype): + for xp in (numpy, cupy): + a = testing.shaped_arange((2,), xp, dtype) + with pytest.raises(xp.linalg.LinAlgError): + xp.linalg.det(a) + + @testing.for_dtypes("fdFD") + def test_det_zero_dim(self, dtype): + for xp in (numpy, cupy): + a = testing.shaped_arange((), xp, dtype) + with pytest.raises(xp.linalg.LinAlgError): + xp.linalg.det(a) + + # TODO: remove skipif when MKLD-16626 is resolved + # _getrf_batch does not raise an error with singular matrices. + # Skip running on cpu because dpnp uses _getrf_batch only on cpu. + @pytest.mark.skipif(is_cpu_device(), reason="MKLD-16626") + @testing.for_dtypes("fdFD") + @testing.numpy_cupy_allclose(rtol=1e-3, atol=1e-4) + def test_det_singular(self, xp, dtype): + a = xp.zeros((2, 3, 3), dtype=dtype) + return xp.linalg.det(a) + + +@pytest.mark.usefixtures("allow_fall_back_on_numpy") +class TestSlogdet(unittest.TestCase): + @testing.for_dtypes("fdFD") + @testing.numpy_cupy_allclose(rtol=1e-3, atol=1e-4) + def test_slogdet(self, xp, dtype): + a = testing.shaped_arange((2, 2), xp, dtype) + 1 + sign, logdet = xp.linalg.slogdet(a) + return sign, logdet + + @testing.for_dtypes("fdFD") + @testing.numpy_cupy_allclose(rtol=1e-3, atol=1e-4) + def test_slogdet_3(self, xp, dtype): + a = testing.shaped_arange((2, 2, 2), xp, dtype) + 1 + sign, logdet = xp.linalg.slogdet(a) + return sign, logdet + + @testing.for_dtypes("fdFD") + @testing.numpy_cupy_allclose(rtol=1e-3, atol=1e-4) + def test_slogdet_4(self, xp, dtype): + a = testing.shaped_arange((2, 2, 2, 2), xp, dtype) + 1 + sign, logdet = xp.linalg.slogdet(a) + return sign, logdet + + @testing.for_dtypes("fdFD") + @testing.numpy_cupy_allclose(rtol=1e-3, atol=1e-4) + def test_slogdet_singular(self, xp, dtype): + a = xp.zeros((3, 3), dtype=dtype) + sign, logdet = xp.linalg.slogdet(a) + return sign, logdet + + @testing.for_dtypes("fdFD") + @testing.numpy_cupy_allclose(rtol=1e-3, atol=1e-4) + def test_slogdet_singular_errstate(self, xp, dtype): + a = xp.zeros((3, 3), dtype=dtype) + # TODO: dpnp has no errstate. Probably to be implemented later + # with cupyx.errstate(linalg="raise"): + # `cupy.linalg.slogdet` internally catches `dev_info < 0` from + # cuSOLVER, which should not affect `dev_info > 0` cases. + sign, logdet = xp.linalg.slogdet(a) + return sign, logdet + + @testing.for_dtypes("fdFD") + def test_slogdet_one_dim(self, dtype): + for xp in (numpy, cupy): + a = testing.shaped_arange((2,), xp, dtype) + with pytest.raises(xp.linalg.LinAlgError): + xp.linalg.slogdet(a)