From 2613795b56cebf3b5ae5c64fca69af9886553247 Mon Sep 17 00:00:00 2001 From: Vladislav Perevezentsev Date: Tue, 26 Sep 2023 12:17:55 +0200 Subject: [PATCH 01/57] Correct return of object type at zero copy --- dpnp/dpnp_container.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/dpnp/dpnp_container.py b/dpnp/dpnp_container.py index faf14c3e97bb..c8c0858df7ea 100644 --- a/dpnp/dpnp_container.py +++ b/dpnp/dpnp_container.py @@ -130,7 +130,7 @@ def asarray( ) # return x1 if dpctl returns a zero copy of x1_obj - if array_obj is x1_obj: + if array_obj is x1_obj and isinstance(x1, dpnp_array): return x1 return dpnp_array(array_obj.shape, buffer=array_obj, order=order) From efdcf2a3268dc83821fa661fa779f45835256f8e Mon Sep 17 00:00:00 2001 From: Vladislav Perevezentsev Date: Tue, 26 Sep 2023 12:19:48 +0200 Subject: [PATCH 02/57] Add tests for gh-1570 --- tests/test_sycl_queue.py | 19 +++++++++++++++++++ tests/test_usm_type.py | 12 ++++++++++++ 2 files changed, 31 insertions(+) diff --git a/tests/test_sycl_queue.py b/tests/test_sycl_queue.py index 48a562cc7989..5711ed7adec0 100644 --- a/tests/test_sycl_queue.py +++ b/tests/test_sycl_queue.py @@ -1,10 +1,12 @@ import dpctl +import dpctl.tensor as dpt import numpy import pytest from dpctl.utils import ExecutionPlacementError from numpy.testing import assert_allclose, assert_array_equal, assert_raises import dpnp +from dpnp.dpnp_array import dpnp_array from .helper import assert_dtype_allclose, get_all_dtypes, is_win_platform @@ -1076,6 +1078,23 @@ def test_array_copy(device, func, device_param, queue_param): assert_sycl_queue_equal(result.sycl_queue, dpnp_data.sycl_queue) +@pytest.mark.parametrize( + "copy", [True, False, None], ids=["True", "False", "None"] +) +@pytest.mark.parametrize( + "device", + valid_devices, + ids=[device.filter_string for device in valid_devices], +) +def test_array_creation_from_dpctl(copy, device): + dpt_data = dpt.ones((3, 3), device=device) + + result = dpnp.array(dpt_data, copy=copy) + + assert_sycl_queue_equal(result.sycl_queue, dpt_data.sycl_queue) + assert isinstance(result, dpnp_array) + + @pytest.mark.parametrize( "device", valid_devices, diff --git a/tests/test_usm_type.py b/tests/test_usm_type.py index 79125a5376b7..bd55670b2a64 100644 --- a/tests/test_usm_type.py +++ b/tests/test_usm_type.py @@ -1,5 +1,6 @@ from math import prod +import dpctl.tensor as dpt import dpctl.utils as du import pytest @@ -180,6 +181,17 @@ def test_array_copy(func, usm_type_x, usm_type_y): assert y.usm_type == usm_type_y +@pytest.mark.parametrize( + "copy", [True, False, None], ids=["True", "False", "None"] +) +@pytest.mark.parametrize("usm_type_x", list_of_usm_types, ids=list_of_usm_types) +def test_array_creation_from_dpctl(copy, usm_type_x): + x = dpt.ones((3, 3), usm_type=usm_type_x) + y = dp.array(x, copy=copy) + + assert y.usm_type == usm_type_x + + @pytest.mark.parametrize( "usm_type_start", list_of_usm_types, ids=list_of_usm_types ) From 6500a6c39006650468db92219ed4cfbf7563d745 Mon Sep 17 00:00:00 2001 From: Vladislav Perevezentsev Date: Thu, 5 Oct 2023 19:34:58 +0200 Subject: [PATCH 03/57] Add dpnp.linalg.solve() function --- dpnp/backend/extensions/lapack/CMakeLists.txt | 1 + dpnp/backend/extensions/lapack/gesv.cpp | 236 ++++++++++++++++++ dpnp/backend/extensions/lapack/gesv.hpp | 50 ++++ dpnp/backend/extensions/lapack/lapack_py.cpp | 10 + .../extensions/lapack/types_matrix.hpp | 26 ++ dpnp/linalg/dpnp_iface_linalg.py | 56 ++++- dpnp/linalg/dpnp_utils_linalg.py | 71 ++++++ 7 files changed, 449 insertions(+), 1 deletion(-) create mode 100644 dpnp/backend/extensions/lapack/gesv.cpp create mode 100644 dpnp/backend/extensions/lapack/gesv.hpp diff --git a/dpnp/backend/extensions/lapack/CMakeLists.txt b/dpnp/backend/extensions/lapack/CMakeLists.txt index 0c90b4f0ca52..81196a78ab98 100644 --- a/dpnp/backend/extensions/lapack/CMakeLists.txt +++ b/dpnp/backend/extensions/lapack/CMakeLists.txt @@ -29,6 +29,7 @@ pybind11_add_module(${python_module_name} MODULE lapack_py.cpp heevd.cpp syevd.cpp + gesv.cpp ) if (WIN32) diff --git a/dpnp/backend/extensions/lapack/gesv.cpp b/dpnp/backend/extensions/lapack/gesv.cpp new file mode 100644 index 000000000000..277886b78954 --- /dev/null +++ b/dpnp/backend/extensions/lapack/gesv.cpp @@ -0,0 +1,236 @@ +//***************************************************************************** +// 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 "gesv.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 (*gesv_impl_fn_ptr_t)(sycl::queue, + const std::int64_t, + const std::int64_t, + char *, + std::int64_t, + std::int64_t *, + char *, + std::int64_t, + std::vector &, + const std::vector &); + +static gesv_impl_fn_ptr_t gesv_dispatch_vector[dpctl_td_ns::num_types]; + +template +static sycl::event gesv_impl(sycl::queue exec_q, + const std::int64_t n, + const std::int64_t nrhs, + char *in_a, + std::int64_t lda, + std::int64_t *ipiv, + char *in_b, + std::int64_t ldb, + std::vector &host_task_events, + const std::vector &depends) +{ + type_utils::validate_type_for_device(exec_q); + + T *a = reinterpret_cast(in_a); + T *b = reinterpret_cast(in_b); + + const std::int64_t scratchpad_size = + mkl_lapack::gesv_scratchpad_size(exec_q, n, nrhs, lda, ldb); + T *scratchpad = nullptr; + + std::stringstream error_msg; + std::int64_t info = 0; + + sycl::event gesv_event; + try { + scratchpad = sycl::malloc_device(scratchpad_size, exec_q); + + gesv_event = mkl_lapack::gesv( + exec_q, + n, // The order of the matrix A (0 ≤ n). + nrhs, // The number of right-hand sides B (0 ≤ nrhs). + a, // Pointer to the square coefficient matrix A (n x n). + lda, // The leading dimension of a, must be at least max(1, n). + ipiv, // The pivot indices that define the permutation matrix P; + // row i of the matrix was interchanged with row ipiv(i), + // must be at least max(1, n). + b, // Pointer to the right hand side matrix B (n x nrhs). + ldb, // The leading dimension of b, must be at least max(1, n). + scratchpad, // Pointer to scratchpad memory to be used by MKL + // routine for storing intermediate results. + scratchpad_size, depends); + } catch (mkl_lapack::exception const &e) { + error_msg + << "Unexpected MKL exception caught during gesv() call:\nreason: " + << e.what() << "\ninfo: " << e.info(); + info = e.info(); + } catch (sycl::exception const &e) { + error_msg << "Unexpected SYCL exception caught during gesv() call:\n" + << e.what(); + info = -1; + } + + if (info != 0) // an unexected error occurs + { + if (scratchpad != nullptr) { + sycl::free(scratchpad, exec_q); + } + throw std::runtime_error(error_msg.str()); + } + + sycl::event clean_up_event = exec_q.submit([&](sycl::handler &cgh) { + cgh.depends_on(gesv_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 gesv_event; +} + +sycl::event gesv(sycl::queue exec_q, + dpctl::tensor::usm_ndarray coeff_matrix, + dpctl::tensor::usm_ndarray hand_sides, + const std::vector &depends) +{ + // check ndim hand_sides + + const py::ssize_t *coeff_matrix_shape = coeff_matrix.get_shape_raw(); + const py::ssize_t *hand_sides_shape = hand_sides.get_shape_raw(); + + if (coeff_matrix_shape[0] != coeff_matrix_shape[1]) { + throw py::value_error("The input coefficients array must be square "); + } + + // elif check shape coeff_matrix and hand_sides + + // check compatibility of execution queue and allocation queue + if (!dpctl::utils::queues_are_compatible(exec_q, + {coeff_matrix, hand_sides})) { + throw py::value_error( + "Execution queue is not compatible with allocation queues"); + } + + auto const &overlap = dpctl::tensor::overlap::MemoryOverlap(); + if (overlap(coeff_matrix, hand_sides)) { + throw py::value_error("Arrays with coefficients and hand sides are " + "overlapping segments of memory"); + } + + bool is_coeff_matrix_f_contig = coeff_matrix.is_f_contiguous(); + bool is_hand_sides_c_contig = hand_sides.is_c_contiguous(); + if (!is_coeff_matrix_f_contig) { + throw py::value_error("An array with coefficients " + "must be F-contiguous"); + } + else if (!is_hand_sides_c_contig) { + throw py::value_error( + "An array with the output solutions of the coefficient matrix" + "must be C-contiguous"); + } + + auto array_types = dpctl_td_ns::usm_ndarray_types(); + int coeff_matrix_type_id = + array_types.typenum_to_lookup_id(coeff_matrix.get_typenum()); + int hand_sides_type_id = + array_types.typenum_to_lookup_id(hand_sides.get_typenum()); + + if (coeff_matrix_type_id != hand_sides_type_id) { + throw py::value_error( + "Types of coefficients and hand sides are missmatched"); + } + + gesv_impl_fn_ptr_t gesv_fn = gesv_dispatch_vector[coeff_matrix_type_id]; + if (gesv_fn == nullptr) { + throw py::value_error("No gesv implementation defined for a type of " + "coefficient matrix and hand sides"); + } + + char *coeff_matrix_data = coeff_matrix.get_data(); + char *hand_sides_data = hand_sides.get_data(); + + const std::int64_t n = coeff_matrix_shape[0]; + const std::int64_t nrhs = hand_sides_shape[0]; + const std::int64_t m = hand_sides_shape[0]; + + const std::int64_t lda = std::max(1UL, n); + const std::int64_t ldb = std::max(1UL, m); + + std::vector ipiv(n); + std::int64_t *d_ipiv = sycl::malloc_device(n, exec_q); + + std::vector host_task_events; + sycl::event gesv_ev = + gesv_fn(exec_q, n, nrhs, coeff_matrix_data, lda, d_ipiv, + hand_sides_data, ldb, host_task_events, depends); + + return gesv_ev; +} + +template +struct GesvContigFactory +{ + fnT get() + { + if constexpr (types::GesvTypePairSupportFactory::is_defined) { + return gesv_impl; + } + else { + return nullptr; + } + } +}; + +void init_gesv_dispatch_vector(void) +{ + dpctl_td_ns::DispatchVectorBuilder + contig; + contig.populate_dispatch_vector(gesv_dispatch_vector); +} +} // namespace lapack +} // namespace ext +} // namespace backend +} // namespace dpnp diff --git a/dpnp/backend/extensions/lapack/gesv.hpp b/dpnp/backend/extensions/lapack/gesv.hpp new file mode 100644 index 000000000000..6c84b08acaf1 --- /dev/null +++ b/dpnp/backend/extensions/lapack/gesv.hpp @@ -0,0 +1,50 @@ +//***************************************************************************** +// 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 sycl::event gesv(sycl::queue exec_q, + dpctl::tensor::usm_ndarray coeff_matrix, + dpctl::tensor::usm_ndarray hand_sides, + const std::vector &depends); + +extern void init_gesv_dispatch_vector(void); +} // 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 97b67d59e24e..e7c419942cdb 100644 --- a/dpnp/backend/extensions/lapack/lapack_py.cpp +++ b/dpnp/backend/extensions/lapack/lapack_py.cpp @@ -30,6 +30,7 @@ #include #include +#include "gesv.hpp" #include "heevd.hpp" #include "syevd.hpp" @@ -40,6 +41,7 @@ namespace py = pybind11; void init_dispatch_vectors(void) { lapack_ext::init_syevd_dispatch_vector(); + lapack_ext::init_gesv_dispatch_vector(); } // populate dispatch tables @@ -66,4 +68,12 @@ PYBIND11_MODULE(_lapack_impl, m) py::arg("sycl_queue"), py::arg("jobz"), py::arg("upper_lower"), py::arg("eig_vecs"), py::arg("eig_vals"), py::arg("depends") = py::list()); + + m.def("_gesv", &lapack_ext::gesv, + "Call `gesv` from OneMKL LAPACK library to return " + "solution to the system of linear equations with a square " + "coefficient matrix A" + "and multiple right-hand sides", + py::arg("sycl_queue"), py::arg("coeff_matrix"), py::arg("hand_sides"), + py::arg("depends") = py::list()); } diff --git a/dpnp/backend/extensions/lapack/types_matrix.hpp b/dpnp/backend/extensions/lapack/types_matrix.hpp index 3cab18d3c63d..9d0a55c41ccf 100644 --- a/dpnp/backend/extensions/lapack/types_matrix.hpp +++ b/dpnp/backend/extensions/lapack/types_matrix.hpp @@ -80,6 +80,32 @@ struct SyevdTypePairSupportFactory // 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::gesv + * function. + * + * @tparam T Type of array containing input matrix A and an output arrays with + * coefficient matrix and hand sides. + */ +template +struct GesvTypePairSupportFactory +{ + 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; +}; } // namespace types } // namespace lapack } // namespace ext diff --git a/dpnp/linalg/dpnp_iface_linalg.py b/dpnp/linalg/dpnp_iface_linalg.py index 3addcdc32585..a089473fffb0 100644 --- a/dpnp/linalg/dpnp_iface_linalg.py +++ b/dpnp/linalg/dpnp_iface_linalg.py @@ -47,7 +47,7 @@ from dpnp.dpnp_utils import * from dpnp.linalg.dpnp_algo_linalg import * -from .dpnp_utils_linalg import dpnp_eigh +from .dpnp_utils_linalg import dpnp_eigh, dpnp_solve __all__ = [ "cholesky", @@ -62,6 +62,7 @@ "multi_dot", "norm", "qr", + "solve", "svd", ] @@ -498,6 +499,59 @@ def qr(x1, mode="reduced"): return call_origin(numpy.linalg.qr, x1, mode) +def solve(a, b): + """ + Solve a linear matrix equation, or system of linear scalar equations. + + For full documentation refer to :obj:`numpy.linalg.solve`. + + Returns + ------- + out : {(…, M,), (…, M, K)} dpnp.ndarray + Solution to the system ax = b. Returned shape is identical to b. + + Limitations + ----------- + Parameter `a` is supported as :class:`dpnp.ndarray` or :class:`dpctl.tensor.usm_ndarray`. + Input array data types are limited by supported DPNP :ref:`Data types`. + + See Also + -------- + :obj:`dpnp.dot` : Returns the dot product of two arrays. + + Examples + -------- + >>> import dpnp as dp + >>> a = dp.array([[1, 2], [3, 5]]) + >>> b = dp.array([1, 2]) + >>> x = dp.linalg.solve(a, b) + >>> x + array([-1., 1.]) + + """ + if not dpnp.is_supported_array_type(a): + raise TypeError( + "An array must be any of supported type, but got {}".format(type(a)) + ) + + if not dpnp.is_supported_array_type(b): + raise TypeError( + "An array must be any of supported type, but got {}".format(type(b)) + ) + + if a.ndim < 2: + raise ValueError( + f"{a.ndim}-dimensional array given. Array must be " + "at least two-dimensional" + ) + + m, n = a.shape[-2:] + if m != n: + raise ValueError("Last 2 dimensions of the array must be square") + + return dpnp_solve(a, b) + + def svd(x1, full_matrices=True, compute_uv=True, hermitian=False): """ Singular Value Decomposition. diff --git a/dpnp/linalg/dpnp_utils_linalg.py b/dpnp/linalg/dpnp_utils_linalg.py index e818835ecbee..baca6d43ce3e 100644 --- a/dpnp/linalg/dpnp_utils_linalg.py +++ b/dpnp/linalg/dpnp_utils_linalg.py @@ -27,6 +27,7 @@ # ***************************************************************************** +import dpctl import dpctl.tensor._tensor_impl as ti import dpnp @@ -164,3 +165,73 @@ def dpnp_eigh(a, UPLO): ht_copy_ev.wait() return w, out_v + + +def dpnp_solve(a, b): + """ + dpnp_solve(a, b) + + Return the the solution to the system of linear equations with + a square coefficient matrix `a` and multiple right-hand sides `b`. + + """ + + a_usm_arr = dpnp.get_usm_ndarray(a) + b_usm_arr = dpnp.get_usm_ndarray(b) + + b_order = "C" if b.flags.c_contiguous else "F" + + if a.dtype != b.dtype: + raise ValueError("a and b must be of the same type") + + exec_q = dpctl.utils.get_execution_queue((a.sycl_queue, b.sycl_queue)) + if exec_q is None: + raise ValueError( + "Execution placement can not be unambiguously inferred " + "from input arguments." + ) + + if dpnp.issubdtype(a.dtype, dpnp.floating): + res_type = ( + a.dtype if exec_q.sycl_device.has_aspect_fp64 else dpnp.float32 + ) + elif dpnp.issubdtype(a.dtype, dpnp.complexfloating): + res_type = ( + a.dtype if exec_q.sycl_device.has_aspect_fp64 else dpnp.complex64 + ) + else: + res_type = ( + dpnp.float64 if exec_q.sycl_device.has_aspect_fp64 else dpnp.float32 + ) + + a_f = dpnp.empty_like(a, order="F", dtype=res_type) + b_f = dpnp.empty_like(b, order="F", dtype=res_type) + + a_ht_copy_ev, a_copy_ev = ti._copy_usm_ndarray_into_usm_ndarray( + src=a_usm_arr, dst=a_f.get_array(), sycl_queue=a.sycl_queue + ) + b_ht_copy_ev, b_copy_ev = ti._copy_usm_ndarray_into_usm_ndarray( + src=b_usm_arr, dst=b_f.get_array(), sycl_queue=b.sycl_queue + ) + + lapack_ev = li._gesv( + exec_q, a_f.get_array(), b_f.get_array(), [a_copy_ev, b_copy_ev] + ) + + if b_order != "F": + out_v = dpnp.empty_like(b_f, order=b_order) + ht_copy_out_ev, _ = ti._copy_usm_ndarray_into_usm_ndarray( + src=b_f.get_array(), + dst=out_v.get_array(), + sycl_queue=b.sycl_queue, + depends=[lapack_ev], + ) + ht_copy_out_ev.wait() + else: + out_v = b_f + + lapack_ev.wait() + b_ht_copy_ev.wait() + a_ht_copy_ev.wait() + + return out_v From 51663d0d038a6c2dac1ba2d79635b88043bcbdfc Mon Sep 17 00:00:00 2001 From: Vladislav Perevezentsev Date: Fri, 6 Oct 2023 13:58:29 +0200 Subject: [PATCH 04/57] Check validity of input array shapes --- dpnp/linalg/dpnp_iface_linalg.py | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/dpnp/linalg/dpnp_iface_linalg.py b/dpnp/linalg/dpnp_iface_linalg.py index a089473fffb0..ae007ee9e823 100644 --- a/dpnp/linalg/dpnp_iface_linalg.py +++ b/dpnp/linalg/dpnp_iface_linalg.py @@ -549,6 +549,15 @@ def solve(a, b): if m != n: raise ValueError("Last 2 dimensions of the array must be square") + if not ( + (a.ndim == b.ndim or a.ndim == b.ndim + 1) + and a.shape[:-1] == b.shape[: a.ndim - 1] + ): + raise ValueError( + "a must have (..., M, M) shape and b must have (..., M) " + "or (..., M, K)" + ) + return dpnp_solve(a, b) From 00b7041e39f024e85b1b46b832d2c916b0f1e44d Mon Sep 17 00:00:00 2001 From: Vladislav Perevezentsev Date: Fri, 6 Oct 2023 14:16:04 +0200 Subject: [PATCH 05/57] Add logic for a.ndim > 2 --- dpnp/linalg/dpnp_utils_linalg.py | 132 ++++++++++++++++++++++++------- 1 file changed, 104 insertions(+), 28 deletions(-) diff --git a/dpnp/linalg/dpnp_utils_linalg.py b/dpnp/linalg/dpnp_utils_linalg.py index baca6d43ce3e..48940c15a193 100644 --- a/dpnp/linalg/dpnp_utils_linalg.py +++ b/dpnp/linalg/dpnp_utils_linalg.py @@ -29,6 +29,7 @@ import dpctl import dpctl.tensor._tensor_impl as ti +from numpy import prod import dpnp import dpnp.backend.extensions.lapack._lapack_impl as li @@ -204,34 +205,109 @@ def dpnp_solve(a, b): dpnp.float64 if exec_q.sycl_device.has_aspect_fp64 else dpnp.float32 ) - a_f = dpnp.empty_like(a, order="F", dtype=res_type) - b_f = dpnp.empty_like(b, order="F", dtype=res_type) - - a_ht_copy_ev, a_copy_ev = ti._copy_usm_ndarray_into_usm_ndarray( - src=a_usm_arr, dst=a_f.get_array(), sycl_queue=a.sycl_queue - ) - b_ht_copy_ev, b_copy_ev = ti._copy_usm_ndarray_into_usm_ndarray( - src=b_usm_arr, dst=b_f.get_array(), sycl_queue=b.sycl_queue - ) - - lapack_ev = li._gesv( - exec_q, a_f.get_array(), b_f.get_array(), [a_copy_ev, b_copy_ev] - ) - - if b_order != "F": - out_v = dpnp.empty_like(b_f, order=b_order) - ht_copy_out_ev, _ = ti._copy_usm_ndarray_into_usm_ndarray( - src=b_f.get_array(), - dst=out_v.get_array(), - sycl_queue=b.sycl_queue, - depends=[lapack_ev], - ) - ht_copy_out_ev.wait() + if a.ndim > 2: + reshape = False + orig_shape_b = b.shape + if a.ndim > 3: + # get 3d input arrays by reshape + if a.ndim == b.ndim: + b = b.reshape(prod(b.shape[:-2]), b.shape[-2], b.shape[-1]) + else: + b = b.reshape(prod(b.shape[:-1]), b.shape[-1]) + + a = a.reshape(prod(a.shape[:-2]), a.shape[-2], a.shape[-1]) + + a_usm_arr = dpnp.get_usm_ndarray(a) + b_usm_arr = dpnp.get_usm_ndarray(b) + reshape = True + + op_count = a.shape[0] + if op_count == 0: + return dpnp.empty_like(b, dtype=res_type) + + coeff_vecs = [None] * op_count + val_vecs = [None] * op_count + a_ht_copy_ev = [None] * op_count + b_ht_copy_ev = [None] * op_count + ht_lapack_ev = [None] * op_count + + for i in range(op_count): + # oneMKL LAPACK assumes fortran-like array as input, so + # allocate a memory with 'F' order for dpnp array of coefficient matrix + # and multiple right-hand sides + coeff_vecs[i] = dpnp.empty_like(a[i], order="F", dtype=res_type) + val_vecs[i] = dpnp.empty_like(b[i], order="F", dtype=res_type) + + # use DPCTL tensor function to fill the array of coefficient matrix + # and multiple right-hand sides with content of input array + a_ht_copy_ev[i], a_copy_ev = ti._copy_usm_ndarray_into_usm_ndarray( + src=a_usm_arr[i], + dst=coeff_vecs[i].get_array(), + sycl_queue=a.sycl_queue, + ) + b_ht_copy_ev[i], b_copy_ev = ti._copy_usm_ndarray_into_usm_ndarray( + src=b_usm_arr[i], + dst=val_vecs[i].get_array(), + sycl_queue=b.sycl_queue, + ) + + # call LAPACK extension function to get the solution of the system of linear + # equations with a portion of the coefficients square matrix + ht_lapack_ev[i] = li._gesv( + exec_q, + coeff_vecs[i].get_array(), + val_vecs[i].get_array(), + depends=[a_copy_ev, b_copy_ev], + ) + + for i in range(op_count): + ht_lapack_ev[i].wait() + b_ht_copy_ev[i].wait() + a_ht_copy_ev[i].wait() + + # combine the list of solutions into a single array + out_v = dpnp.array(val_vecs, order=b_order) + if reshape: + # shape of the out_t must be equal to the shape of the right-hand sides + out_v = out_v.reshape(orig_shape_b) + return out_v else: - out_v = b_f + # oneMKL LAPACK assumes fortran-like array as input, so + # allocate a memory with 'F' order for dpnp array of coefficient matrix + # and multiple right-hand sides + a_f = dpnp.empty_like(a, order="F", dtype=res_type) + b_f = dpnp.empty_like(b, order="F", dtype=res_type) + + # use DPCTL tensor function to fill the array of coefficient matrix + # and multiple right-hand sides with content of input array + a_ht_copy_ev, a_copy_ev = ti._copy_usm_ndarray_into_usm_ndarray( + src=a_usm_arr, dst=a_f.get_array(), sycl_queue=a.sycl_queue + ) + b_ht_copy_ev, b_copy_ev = ti._copy_usm_ndarray_into_usm_ndarray( + src=b_usm_arr, dst=b_f.get_array(), sycl_queue=b.sycl_queue + ) + + # call LAPACK extension function to get the solution of the system of linear + # equations with the coefficients square matrix + lapack_ev = li._gesv( + exec_q, a_f.get_array(), b_f.get_array(), [a_copy_ev, b_copy_ev] + ) + + if b_order != "F": + # need to align order of the result of solutions with the right-hand sides + out_v = dpnp.empty_like(b_f, order=b_order) + ht_copy_out_ev, _ = ti._copy_usm_ndarray_into_usm_ndarray( + src=b_f.get_array(), + dst=out_v.get_array(), + sycl_queue=b.sycl_queue, + depends=[lapack_ev], + ) + ht_copy_out_ev.wait() + else: + out_v = b_f - lapack_ev.wait() - b_ht_copy_ev.wait() - a_ht_copy_ev.wait() + lapack_ev.wait() + b_ht_copy_ev.wait() + a_ht_copy_ev.wait() - return out_v + return out_v From e4e15ad2d2faad28cf5169ce5bda21b9c7118e7c Mon Sep 17 00:00:00 2001 From: Vladislav Perevezentsev Date: Fri, 6 Oct 2023 14:19:50 +0200 Subject: [PATCH 06/57] Raise value_error if coeff_matrix_nd != 2 in gesv --- dpnp/backend/extensions/lapack/gesv.cpp | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/dpnp/backend/extensions/lapack/gesv.cpp b/dpnp/backend/extensions/lapack/gesv.cpp index 277886b78954..d84f51a0005e 100644 --- a/dpnp/backend/extensions/lapack/gesv.cpp +++ b/dpnp/backend/extensions/lapack/gesv.cpp @@ -135,7 +135,13 @@ sycl::event gesv(sycl::queue exec_q, dpctl::tensor::usm_ndarray hand_sides, const std::vector &depends) { - // check ndim hand_sides + const int coeff_matrix_nd = coeff_matrix.get_ndim(); + + if (coeff_matrix_nd != 2) { + throw py::value_error( + "Unexpected ndim=" + std::to_string(coeff_matrix_nd) + + " of an input array with coefficients"); + } const py::ssize_t *coeff_matrix_shape = coeff_matrix.get_shape_raw(); const py::ssize_t *hand_sides_shape = hand_sides.get_shape_raw(); @@ -160,16 +166,10 @@ sycl::event gesv(sycl::queue exec_q, } bool is_coeff_matrix_f_contig = coeff_matrix.is_f_contiguous(); - bool is_hand_sides_c_contig = hand_sides.is_c_contiguous(); if (!is_coeff_matrix_f_contig) { throw py::value_error("An array with coefficients " "must be F-contiguous"); } - else if (!is_hand_sides_c_contig) { - throw py::value_error( - "An array with the output solutions of the coefficient matrix" - "must be C-contiguous"); - } auto array_types = dpctl_td_ns::usm_ndarray_types(); int coeff_matrix_type_id = From e762c66fa0fecd718868d8df388f0274f64c63fc Mon Sep 17 00:00:00 2001 From: Vladislav Perevezentsev Date: Fri, 6 Oct 2023 14:28:38 +0200 Subject: [PATCH 07/57] Add cupy tests for dpnp.linalg.solve() --- .../cupy/linalg_tests/test_solve.py | 80 +++++++++++++++++++ 1 file changed, 80 insertions(+) create mode 100644 tests/third_party/cupy/linalg_tests/test_solve.py diff --git a/tests/third_party/cupy/linalg_tests/test_solve.py b/tests/third_party/cupy/linalg_tests/test_solve.py new file mode 100644 index 000000000000..17cc11be461e --- /dev/null +++ b/tests/third_party/cupy/linalg_tests/test_solve.py @@ -0,0 +1,80 @@ +import unittest + +import numpy +import pytest + +import dpnp as cupy +from tests.helper import has_support_aspect64 +from tests.third_party.cupy import testing + + +@testing.parameterize( + *testing.product( + { + "order": ["C", "F"], + } + ) +) +class TestSolve(unittest.TestCase): + # TODO: add get_batched_gesv_limit + # def setUp(self): + # if self.batched_gesv_limit is not None: + # self.old_limit = get_batched_gesv_limit() + # set_batched_gesv_limit(self.batched_gesv_limit) + + # def tearDown(self): + # if self.batched_gesv_limit is not None: + # set_batched_gesv_limit(self.old_limit) + + @testing.for_dtypes("ifdFD") + @testing.numpy_cupy_allclose( + atol=1e-3, contiguous_check=False, type_check=has_support_aspect64() + ) + def check_x(self, a_shape, b_shape, xp, dtype): + a = testing.shaped_random(a_shape, xp, dtype=dtype, seed=0, scale=20) + b = testing.shaped_random(b_shape, xp, dtype=dtype, seed=1) + a = a.copy(order=self.order) + b = b.copy(order=self.order) + a_copy = a.copy() + b_copy = b.copy() + result = xp.linalg.solve(a, b) + numpy.testing.assert_array_equal(a_copy, a) + numpy.testing.assert_array_equal(b_copy, b) + return result + + def test_solve(self): + self.check_x((4, 4), (4,)) + self.check_x((5, 5), (5, 2)) + self.check_x((2, 4, 4), (2, 4)) + self.check_x((2, 5, 5), (2, 5, 2)) + self.check_x((2, 3, 2, 2), (2, 3, 2)) + self.check_x((2, 3, 3, 3), (2, 3, 3, 2)) + self.check_x((0, 0), (0,)) + self.check_x((0, 0), (0, 2)) + self.check_x((0, 2, 2), (0, 2)) + self.check_x((0, 2, 2), (0, 2, 3)) + + def check_shape(self, a_shape, b_shape, error_type): + for xp in (numpy, cupy): + a = xp.random.rand(*a_shape) + b = xp.random.rand(*b_shape) + with pytest.raises(error_type): + xp.linalg.solve(a, b) + + # dpnp.linalg.solve() raises RuntimeError instead of numpy.linalg.LinAlgError + # @testing.numpy_cupy_allclose() + # def test_solve_singular_empty(self, xp): + # a = xp.zeros((3, 3)) # singular + # b = xp.empty((3, 0)) # nrhs = 0 + # # LinAlgError("Singular matrix") is not raised + # return xp.linalg.solve(a, b) + + # dpnp.linalg.solve() raises RuntimeError instead of numpy.linalg.LinAlgError + def test_invalid_shape(self): + # self.check_shape((2, 3), (4,), numpy.linalg.LinAlgError) + self.check_shape((3, 3), (2,), ValueError) + self.check_shape((3, 3), (2, 2), ValueError) + # self.check_shape((3, 3, 4), (3,), numpy.linalg.LinAlgError) + self.check_shape((2, 3, 3), (3,), ValueError) + self.check_shape((3, 3), (0,), ValueError) + # self.check_shape((0, 3, 4), (3,), numpy.linalg.LinAlgError) From a0f76d529d12c57aea5fddeae0066bb655edb08d Mon Sep 17 00:00:00 2001 From: Vladislav Perevezentsev Date: Tue, 10 Oct 2023 22:34:46 +0200 Subject: [PATCH 08/57] Add LinAlgError exception and extend error handling for mkl::lapack::gesv --- dpnp/backend/extensions/lapack/gesv.cpp | 46 +++++++++++++++++-- dpnp/backend/extensions/lapack/lapack_py.cpp | 4 ++ .../extensions/lapack/linalg_exceptions.hpp | 29 ++++++++++++ 3 files changed, 75 insertions(+), 4 deletions(-) create mode 100644 dpnp/backend/extensions/lapack/linalg_exceptions.hpp diff --git a/dpnp/backend/extensions/lapack/gesv.cpp b/dpnp/backend/extensions/lapack/gesv.cpp index d84f51a0005e..87a9941df879 100644 --- a/dpnp/backend/extensions/lapack/gesv.cpp +++ b/dpnp/backend/extensions/lapack/gesv.cpp @@ -30,6 +30,7 @@ #include "utils/type_utils.hpp" #include "gesv.hpp" +#include "linalg_exceptions.hpp" #include "types_matrix.hpp" #include "dpnp_utils.hpp" @@ -102,14 +103,51 @@ static sycl::event gesv_impl(sycl::queue exec_q, // routine for storing intermediate results. scratchpad_size, depends); } catch (mkl_lapack::exception const &e) { - error_msg - << "Unexpected MKL exception caught during gesv() call:\nreason: " - << e.what() << "\ninfo: " << e.info(); info = e.info(); + + if (info < 0) { + error_msg << "Parameter number " << -info + << " had an illegal value."; + } + else if (info > 0) { + T host_U; + exec_q.memcpy(&host_U, &a[(info - 1) * lda + info - 1], sizeof(T)) + .wait(); + + using ThresholdType = typename std::conditional< + std::is_same::value, float, + typename std::conditional< + std::is_same::value, double, + typename std::conditional< + std::is_same>::value, float, + double>::type>::type>::type; + + const auto threshold = + std::numeric_limits::epsilon() * 100; + if (std::abs(host_U) < threshold) { + sycl::free(scratchpad, exec_q); + throw LinAlgError("The input coefficient matrix is singular."); + } + else { + error_msg << "Unexpected MKL exception caught during gesv() " + "call:\nreason: " + << 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) { error_msg << "Unexpected SYCL exception caught during gesv() call:\n" << e.what(); - info = -1; + info = -11; } if (info != 0) // an unexected error occurs diff --git a/dpnp/backend/extensions/lapack/lapack_py.cpp b/dpnp/backend/extensions/lapack/lapack_py.cpp index e7c419942cdb..8a9c07c20c12 100644 --- a/dpnp/backend/extensions/lapack/lapack_py.cpp +++ b/dpnp/backend/extensions/lapack/lapack_py.cpp @@ -32,6 +32,7 @@ #include "gesv.hpp" #include "heevd.hpp" +#include "linalg_exceptions.hpp" #include "syevd.hpp" namespace lapack_ext = dpnp::backend::ext::lapack; @@ -52,6 +53,9 @@ void init_dispatch_tables(void) PYBIND11_MODULE(_lapack_impl, m) { + py::register_local_exception(m, "LinAlgError", + PyExc_ValueError); + init_dispatch_vectors(); init_dispatch_tables(); diff --git a/dpnp/backend/extensions/lapack/linalg_exceptions.hpp b/dpnp/backend/extensions/lapack/linalg_exceptions.hpp new file mode 100644 index 000000000000..9963b2aa834e --- /dev/null +++ b/dpnp/backend/extensions/lapack/linalg_exceptions.hpp @@ -0,0 +1,29 @@ +#pragma once +#include +#include + +namespace dpnp +{ +namespace backend +{ +namespace ext +{ +namespace lapack +{ +class LinAlgError : public std::exception +{ +public: + explicit LinAlgError(const char *message) : msg_(message) {} + + const char *what() const noexcept override + { + return msg_.c_str(); + } + +private: + std::string msg_; +}; +} // namespace lapack +} // namespace ext +} // namespace backend +} // namespace dpnp From 89f47c75c42c239963e4431a72199b2a8b3fda97 Mon Sep 17 00:00:00 2001 From: Vladislav Perevezentsev Date: Tue, 10 Oct 2023 22:36:00 +0200 Subject: [PATCH 09/57] Update test_solve --- .../cupy/linalg_tests/test_solve.py | 26 +++++++++++-------- 1 file changed, 15 insertions(+), 11 deletions(-) diff --git a/tests/third_party/cupy/linalg_tests/test_solve.py b/tests/third_party/cupy/linalg_tests/test_solve.py index 17cc11be461e..68c4a0ea1b60 100644 --- a/tests/third_party/cupy/linalg_tests/test_solve.py +++ b/tests/third_party/cupy/linalg_tests/test_solve.py @@ -61,20 +61,24 @@ def check_shape(self, a_shape, b_shape, error_type): with pytest.raises(error_type): xp.linalg.solve(a, b) - # dpnp.linalg.solve() raises RuntimeError instead of numpy.linalg.LinAlgError - # @testing.numpy_cupy_allclose() - # def test_solve_singular_empty(self, xp): - # a = xp.zeros((3, 3)) # singular - # b = xp.empty((3, 0)) # nrhs = 0 - # # LinAlgError("Singular matrix") is not raised - # return xp.linalg.solve(a, b) + def test_solve_singular_empty(self): + for xp in (numpy, cupy): + a = xp.zeros((3, 3)) # singular + b = xp.empty((3, 0)) # nrhs = 0 + with pytest.raises((numpy.linalg.LinAlgError, ValueError)): + xp.linalg.solve(a, b) - # dpnp.linalg.solve() raises RuntimeError instead of numpy.linalg.LinAlgError + # dpnp.linalg.solve() raises a LinAlgError which is defined + # through a ValueError in the C++ bindings using pybind11 def test_invalid_shape(self): - # self.check_shape((2, 3), (4,), numpy.linalg.LinAlgError) + self.check_shape((2, 3), (4,), (numpy.linalg.LinAlgError, ValueError)) self.check_shape((3, 3), (2,), ValueError) self.check_shape((3, 3), (2, 2), ValueError) - # self.check_shape((3, 3, 4), (3,), numpy.linalg.LinAlgError) + self.check_shape( + (3, 3, 4), (3,), (numpy.linalg.LinAlgError, ValueError) + ) self.check_shape((2, 3, 3), (3,), ValueError) self.check_shape((3, 3), (0,), ValueError) - # self.check_shape((0, 3, 4), (3,), numpy.linalg.LinAlgError) + self.check_shape( + (0, 3, 4), (3,), (numpy.linalg.LinAlgError, ValueError) + ) From 159c46019189de5900b8abb7d8344bc239dd7ac6 Mon Sep 17 00:00:00 2001 From: Vladislav Perevezentsev Date: Tue, 10 Oct 2023 22:45:24 +0200 Subject: [PATCH 10/57] Add test_solve to test scope --- .github/workflows/conda-package.yml | 1 + 1 file changed, 1 insertion(+) diff --git a/.github/workflows/conda-package.yml b/.github/workflows/conda-package.yml index 365d8490ec67..f0c28edaaa58 100644 --- a/.github/workflows/conda-package.yml +++ b/.github/workflows/conda-package.yml @@ -29,6 +29,7 @@ env: test_usm_type.py third_party/cupy/core_tests/test_ndarray_complex_ops.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 third_party/cupy/logic_tests/test_truth.py third_party/cupy/manipulation_tests/test_basic.py From 2ad3aa816867f9ef842a50e757d6c1306bdf7534 Mon Sep 17 00:00:00 2001 From: Vladislav Perevezentsev Date: Thu, 12 Oct 2023 16:39:06 +0200 Subject: [PATCH 11/57] Fix getting nrhs to avoid CPU falling tests --- dpnp/backend/extensions/lapack/gesv.cpp | 5 ++--- dpnp/linalg/dpnp_utils_linalg.py | 3 +++ tests/third_party/cupy/linalg_tests/test_solve.py | 10 ++++++++-- 3 files changed, 13 insertions(+), 5 deletions(-) diff --git a/dpnp/backend/extensions/lapack/gesv.cpp b/dpnp/backend/extensions/lapack/gesv.cpp index 87a9941df879..5b8b781d7d4b 100644 --- a/dpnp/backend/extensions/lapack/gesv.cpp +++ b/dpnp/backend/extensions/lapack/gesv.cpp @@ -174,6 +174,7 @@ sycl::event gesv(sycl::queue exec_q, const std::vector &depends) { const int coeff_matrix_nd = coeff_matrix.get_ndim(); + const int hand_sides_nd = hand_sides.get_ndim(); if (coeff_matrix_nd != 2) { throw py::value_error( @@ -188,8 +189,6 @@ sycl::event gesv(sycl::queue exec_q, throw py::value_error("The input coefficients array must be square "); } - // elif check shape coeff_matrix and hand_sides - // check compatibility of execution queue and allocation queue if (!dpctl::utils::queues_are_compatible(exec_q, {coeff_matrix, hand_sides})) { @@ -230,8 +229,8 @@ sycl::event gesv(sycl::queue exec_q, char *hand_sides_data = hand_sides.get_data(); const std::int64_t n = coeff_matrix_shape[0]; - const std::int64_t nrhs = hand_sides_shape[0]; const std::int64_t m = hand_sides_shape[0]; + const std::int64_t nrhs = (hand_sides_nd > 1) ? hand_sides_shape[1] : 1; const std::int64_t lda = std::max(1UL, n); const std::int64_t ldb = std::max(1UL, m); diff --git a/dpnp/linalg/dpnp_utils_linalg.py b/dpnp/linalg/dpnp_utils_linalg.py index 48940c15a193..d6bcd92fb788 100644 --- a/dpnp/linalg/dpnp_utils_linalg.py +++ b/dpnp/linalg/dpnp_utils_linalg.py @@ -205,6 +205,9 @@ def dpnp_solve(a, b): dpnp.float64 if exec_q.sycl_device.has_aspect_fp64 else dpnp.float32 ) + if b.size == 0: + return dpnp.empty(b.shape, dtype=res_type) + if a.ndim > 2: reshape = False orig_shape_b = b.shape diff --git a/tests/third_party/cupy/linalg_tests/test_solve.py b/tests/third_party/cupy/linalg_tests/test_solve.py index 68c4a0ea1b60..6ebb1f213b69 100644 --- a/tests/third_party/cupy/linalg_tests/test_solve.py +++ b/tests/third_party/cupy/linalg_tests/test_solve.py @@ -65,8 +65,14 @@ def test_solve_singular_empty(self): for xp in (numpy, cupy): a = xp.zeros((3, 3)) # singular b = xp.empty((3, 0)) # nrhs = 0 - with pytest.raises((numpy.linalg.LinAlgError, ValueError)): - xp.linalg.solve(a, b) + # numpy <= 1.24.* raises LinAlgError when b.size == 0 + # numpy >= 1.25 returns an empty array + if xp == numpy: + with pytest.raises(numpy.linalg.LinAlgError): + xp.linalg.solve(a, b) + else: + result = xp.linalg.solve(a, b) + assert result.size == 0 # dpnp.linalg.solve() raises a LinAlgError which is defined # through a ValueError in the C++ bindings using pybind11 From ac02cf2368b43da5f5e800bc69252d32abe12140 Mon Sep 17 00:00:00 2001 From: Vladislav Perevezentsev Date: Thu, 12 Oct 2023 17:00:44 +0200 Subject: [PATCH 12/57] Add test_solve to test_sycl_queue --- tests/test_sycl_queue.py | 24 ++++++++++++++++++++++++ 1 file changed, 24 insertions(+) diff --git a/tests/test_sycl_queue.py b/tests/test_sycl_queue.py index a11880d9d444..b8d94abb54b9 100644 --- a/tests/test_sycl_queue.py +++ b/tests/test_sycl_queue.py @@ -1194,3 +1194,27 @@ def test_take(device): result_queue = result.get_array().sycl_queue 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_solve(device): + x = [[1.0, 2.0], [3.0, 5.0]] + y = [1.0, 2.0] + + numpy_x = numpy.array(x) + numpy_y = numpy.array(y) + dpnp_x = dpnp.array(x, device=device) + dpnp_y = dpnp.array(y, device=device) + + result = dpnp.linalg.solve(dpnp_x, dpnp_y) + expected = numpy.linalg.solve(numpy_x, numpy_y) + assert_allclose(expected, result, rtol=1e-06) + + result_queue = result.sycl_queue + + assert_sycl_queue_equal(result_queue, dpnp_x.sycl_queue) + assert_sycl_queue_equal(result_queue, dpnp_y.sycl_queue) From 035d9831d3a31898d15242dd11fd6d0db03a5d63 Mon Sep 17 00:00:00 2001 From: Vladislav Perevezentsev Date: Fri, 13 Oct 2023 00:15:41 +0200 Subject: [PATCH 13/57] Add more tests for solve() --- tests/test_linalg.py | 39 ++++++++++++++++++- .../cupy/linalg_tests/test_solve.py | 2 +- 2 files changed, 39 insertions(+), 2 deletions(-) diff --git a/tests/test_linalg.py b/tests/test_linalg.py index 2327ca2e9401..2e9ea5986d8c 100644 --- a/tests/test_linalg.py +++ b/tests/test_linalg.py @@ -1,7 +1,7 @@ import dpctl import numpy import pytest -from numpy.testing import assert_allclose, assert_array_equal +from numpy.testing import assert_allclose, assert_array_equal, assert_raises import dpnp as inp @@ -446,3 +446,40 @@ def test_svd(type, shape): assert_allclose( inp.asnumpy(dpnp_vt)[i, :], np_vt[i, :], rtol=tol, atol=tol ) + + +class TestSolve: + @pytest.mark.parametrize("dtype", get_all_dtypes(no_bool=True)) + def test_solve(self, dtype): + a_np = numpy.array([[1, 0.5], [0.5, 1]]) + a_dp = inp.array(a_np, dtype=dtype) + + expected = numpy.linalg.solve(a_np, a_np) + result = inp.linalg.solve(a_dp, a_dp) + + assert_allclose(expected, result, rtol=1e-06) + + def test_solve_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], + ] + ) + b_np = numpy.array([5, 8, 9, 2, 1]) + + a_dp = inp.array(a_np) + b_dp = inp.array(b_np) + + # positive strides + expected = numpy.linalg.solve(a_np[::2, ::2], b_np[::2]) + result = inp.linalg.solve(a_dp[::2, ::2], b_dp[::2]) + assert_allclose(expected, result, rtol=1e-06) + + # negative strides + expected = numpy.linalg.solve(a_np[::-2, ::-2], b_np[::-2]) + result = inp.linalg.solve(a_dp[::-2, ::-2], b_dp[::-2]) + assert_allclose(expected, result, rtol=1e-06) diff --git a/tests/third_party/cupy/linalg_tests/test_solve.py b/tests/third_party/cupy/linalg_tests/test_solve.py index 6ebb1f213b69..2af6362377ca 100644 --- a/tests/third_party/cupy/linalg_tests/test_solve.py +++ b/tests/third_party/cupy/linalg_tests/test_solve.py @@ -67,7 +67,7 @@ def test_solve_singular_empty(self): b = xp.empty((3, 0)) # nrhs = 0 # numpy <= 1.24.* raises LinAlgError when b.size == 0 # numpy >= 1.25 returns an empty array - if xp == numpy: + if xp == numpy and testing.numpy_satisfies("<1.25.0"): with pytest.raises(numpy.linalg.LinAlgError): xp.linalg.solve(a, b) else: From f80ba4b6844d6a6672082c9923766b3e051d5edd Mon Sep 17 00:00:00 2001 From: Vladislav Perevezentsev Date: Fri, 13 Oct 2023 14:01:07 +0200 Subject: [PATCH 14/57] Register a LinAlgError in dpnp.linalg submodule --- dpnp/backend/extensions/lapack/lapack_py.cpp | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/dpnp/backend/extensions/lapack/lapack_py.cpp b/dpnp/backend/extensions/lapack/lapack_py.cpp index 8a9c07c20c12..5ae26d595cb8 100644 --- a/dpnp/backend/extensions/lapack/lapack_py.cpp +++ b/dpnp/backend/extensions/lapack/lapack_py.cpp @@ -53,8 +53,9 @@ void init_dispatch_tables(void) PYBIND11_MODULE(_lapack_impl, m) { - py::register_local_exception(m, "LinAlgError", - PyExc_ValueError); + py::module_ linalg_module = py::module_::import("dpnp.linalg"); + py::register_exception( + linalg_module, "LinAlgError", PyExc_ValueError); init_dispatch_vectors(); init_dispatch_tables(); From 9a3db74320da194594f8fc69dea1d99f85ff9708 Mon Sep 17 00:00:00 2001 From: Vladislav Perevezentsev Date: Fri, 13 Oct 2023 15:11:04 +0200 Subject: [PATCH 15/57] Raise dpnp.linalg.LinAlgError in solve() --- dpnp/linalg/dpnp_iface_linalg.py | 8 +++-- tests/test_linalg.py | 4 +-- .../cupy/linalg_tests/test_solve.py | 33 ++++++++++--------- 3 files changed, 25 insertions(+), 20 deletions(-) diff --git a/dpnp/linalg/dpnp_iface_linalg.py b/dpnp/linalg/dpnp_iface_linalg.py index ae007ee9e823..b9a0d889e046 100644 --- a/dpnp/linalg/dpnp_iface_linalg.py +++ b/dpnp/linalg/dpnp_iface_linalg.py @@ -540,20 +540,22 @@ def solve(a, b): ) if a.ndim < 2: - raise ValueError( + raise dpnp.linalg.LinAlgError( f"{a.ndim}-dimensional array given. Array must be " "at least two-dimensional" ) m, n = a.shape[-2:] if m != n: - raise ValueError("Last 2 dimensions of the array must be square") + raise dpnp.linalg.LinAlgError( + "Last 2 dimensions of the array must be square" + ) if not ( (a.ndim == b.ndim or a.ndim == b.ndim + 1) and a.shape[:-1] == b.shape[: a.ndim - 1] ): - raise ValueError( + raise dpnp.linalg.LinAlgError( "a must have (..., M, M) shape and b must have (..., M) " "or (..., M, K)" ) diff --git a/tests/test_linalg.py b/tests/test_linalg.py index 2e9ea5986d8c..7782fde2431c 100644 --- a/tests/test_linalg.py +++ b/tests/test_linalg.py @@ -477,9 +477,9 @@ def test_solve_strides(self): # positive strides expected = numpy.linalg.solve(a_np[::2, ::2], b_np[::2]) result = inp.linalg.solve(a_dp[::2, ::2], b_dp[::2]) - assert_allclose(expected, result, rtol=1e-06) + assert_allclose(expected, result, rtol=1e-05) # negative strides expected = numpy.linalg.solve(a_np[::-2, ::-2], b_np[::-2]) result = inp.linalg.solve(a_dp[::-2, ::-2], b_dp[::-2]) - assert_allclose(expected, result, rtol=1e-06) + assert_allclose(expected, result, rtol=1e-05) diff --git a/tests/third_party/cupy/linalg_tests/test_solve.py b/tests/third_party/cupy/linalg_tests/test_solve.py index 2af6362377ca..5de4aacc789b 100644 --- a/tests/third_party/cupy/linalg_tests/test_solve.py +++ b/tests/third_party/cupy/linalg_tests/test_solve.py @@ -54,8 +54,8 @@ def test_solve(self): self.check_x((0, 2, 2), (0, 2)) self.check_x((0, 2, 2), (0, 2, 3)) - def check_shape(self, a_shape, b_shape, error_type): - for xp in (numpy, cupy): + def check_shape(self, a_shape, b_shape, error_types): + for xp, error_type in error_types.items(): a = xp.random.rand(*a_shape) b = xp.random.rand(*b_shape) with pytest.raises(error_type): @@ -74,17 +74,20 @@ def test_solve_singular_empty(self): result = xp.linalg.solve(a, b) assert result.size == 0 - # dpnp.linalg.solve() raises a LinAlgError which is defined - # through a ValueError in the C++ bindings using pybind11 def test_invalid_shape(self): - self.check_shape((2, 3), (4,), (numpy.linalg.LinAlgError, ValueError)) - self.check_shape((3, 3), (2,), ValueError) - self.check_shape((3, 3), (2, 2), ValueError) - self.check_shape( - (3, 3, 4), (3,), (numpy.linalg.LinAlgError, ValueError) - ) - self.check_shape((2, 3, 3), (3,), ValueError) - self.check_shape((3, 3), (0,), ValueError) - self.check_shape( - (0, 3, 4), (3,), (numpy.linalg.LinAlgError, ValueError) - ) + linalg_errors = { + numpy: numpy.linalg.LinAlgError, + cupy: cupy.linalg.LinAlgError, + } + value_errors = { + numpy: ValueError, + cupy: ValueError, + } + + self.check_shape((2, 3), (4,), linalg_errors) + self.check_shape((3, 3), (2,), value_errors) + self.check_shape((3, 3), (2, 2), value_errors) + self.check_shape((3, 3, 4), (3,), linalg_errors) + self.check_shape((2, 3, 3), (3,), value_errors) + self.check_shape((3, 3), (0,), value_errors) + self.check_shape((0, 3, 4), (3,), linalg_errors) From a8b4fec4d8d295bd2c1d8aecfef2f9a8bd51f1f4 Mon Sep 17 00:00:00 2001 From: Vladislav Perevezentsev Date: Sun, 15 Oct 2023 20:01:50 +0200 Subject: [PATCH 16/57] Small changes to the docstrings --- dpnp/backend/extensions/lapack/gesv.cpp | 53 +++++++++++-------- dpnp/backend/extensions/lapack/gesv.hpp | 2 +- dpnp/backend/extensions/lapack/lapack_py.cpp | 10 ++-- .../extensions/lapack/linalg_exceptions.hpp | 25 +++++++++ dpnp/linalg/dpnp_iface_linalg.py | 6 +-- dpnp/linalg/dpnp_utils_linalg.py | 32 ++++++----- 6 files changed, 83 insertions(+), 45 deletions(-) diff --git a/dpnp/backend/extensions/lapack/gesv.cpp b/dpnp/backend/extensions/lapack/gesv.cpp index 5b8b781d7d4b..94b4f326eee1 100644 --- a/dpnp/backend/extensions/lapack/gesv.cpp +++ b/dpnp/backend/extensions/lapack/gesv.cpp @@ -170,67 +170,74 @@ static sycl::event gesv_impl(sycl::queue exec_q, sycl::event gesv(sycl::queue exec_q, dpctl::tensor::usm_ndarray coeff_matrix, - dpctl::tensor::usm_ndarray hand_sides, + dpctl::tensor::usm_ndarray dependent_vals, const std::vector &depends) { const int coeff_matrix_nd = coeff_matrix.get_ndim(); - const int hand_sides_nd = hand_sides.get_ndim(); + const int dependent_vals_nd = dependent_vals.get_ndim(); if (coeff_matrix_nd != 2) { - throw py::value_error( - "Unexpected ndim=" + std::to_string(coeff_matrix_nd) + - " of an input array with coefficients"); + throw py::value_error("The coefficient matrix has ndim=" + + std::to_string(coeff_matrix_nd) + + ", but a 2-dimensional array is expected."); } const py::ssize_t *coeff_matrix_shape = coeff_matrix.get_shape_raw(); - const py::ssize_t *hand_sides_shape = hand_sides.get_shape_raw(); + const py::ssize_t *dependent_vals_shape = dependent_vals.get_shape_raw(); if (coeff_matrix_shape[0] != coeff_matrix_shape[1]) { - throw py::value_error("The input coefficients array must be square "); + throw py::value_error("The coefficient matrix must be square," + " but got a shape of (" + + std::to_string(coeff_matrix_shape[0]) + ", " + + std::to_string(coeff_matrix_shape[1]) + ")."); } // check compatibility of execution queue and allocation queue if (!dpctl::utils::queues_are_compatible(exec_q, - {coeff_matrix, hand_sides})) { + {coeff_matrix, dependent_vals})) + { throw py::value_error( "Execution queue is not compatible with allocation queues"); } auto const &overlap = dpctl::tensor::overlap::MemoryOverlap(); - if (overlap(coeff_matrix, hand_sides)) { - throw py::value_error("Arrays with coefficients and hand sides are " - "overlapping segments of memory"); + if (overlap(coeff_matrix, dependent_vals)) { + throw py::value_error( + "The arrays of coefficients and dependent variables " + "are overlapping segments of memory"); } bool is_coeff_matrix_f_contig = coeff_matrix.is_f_contiguous(); if (!is_coeff_matrix_f_contig) { - throw py::value_error("An array with coefficients " + throw py::value_error("The coefficient matrix " "must be F-contiguous"); } auto array_types = dpctl_td_ns::usm_ndarray_types(); int coeff_matrix_type_id = array_types.typenum_to_lookup_id(coeff_matrix.get_typenum()); - int hand_sides_type_id = - array_types.typenum_to_lookup_id(hand_sides.get_typenum()); + int dependent_vals_type_id = + array_types.typenum_to_lookup_id(dependent_vals.get_typenum()); - if (coeff_matrix_type_id != hand_sides_type_id) { - throw py::value_error( - "Types of coefficients and hand sides are missmatched"); + if (coeff_matrix_type_id != dependent_vals_type_id) { + throw py::value_error("The types of the coefficient matrix and " + "dependent variables are mismatched"); } gesv_impl_fn_ptr_t gesv_fn = gesv_dispatch_vector[coeff_matrix_type_id]; if (gesv_fn == nullptr) { - throw py::value_error("No gesv implementation defined for a type of " - "coefficient matrix and hand sides"); + throw py::value_error( + "No gesv implementation defined for the provided type " + "of the coefficient matrix."); } char *coeff_matrix_data = coeff_matrix.get_data(); - char *hand_sides_data = hand_sides.get_data(); + char *dependent_vals_data = dependent_vals.get_data(); const std::int64_t n = coeff_matrix_shape[0]; - const std::int64_t m = hand_sides_shape[0]; - const std::int64_t nrhs = (hand_sides_nd > 1) ? hand_sides_shape[1] : 1; + const std::int64_t m = dependent_vals_shape[0]; + const std::int64_t nrhs = + (dependent_vals_nd > 1) ? dependent_vals_shape[1] : 1; const std::int64_t lda = std::max(1UL, n); const std::int64_t ldb = std::max(1UL, m); @@ -241,7 +248,7 @@ sycl::event gesv(sycl::queue exec_q, std::vector host_task_events; sycl::event gesv_ev = gesv_fn(exec_q, n, nrhs, coeff_matrix_data, lda, d_ipiv, - hand_sides_data, ldb, host_task_events, depends); + dependent_vals_data, ldb, host_task_events, depends); return gesv_ev; } diff --git a/dpnp/backend/extensions/lapack/gesv.hpp b/dpnp/backend/extensions/lapack/gesv.hpp index 6c84b08acaf1..3b9a099e9795 100644 --- a/dpnp/backend/extensions/lapack/gesv.hpp +++ b/dpnp/backend/extensions/lapack/gesv.hpp @@ -40,7 +40,7 @@ namespace lapack { extern sycl::event gesv(sycl::queue exec_q, dpctl::tensor::usm_ndarray coeff_matrix, - dpctl::tensor::usm_ndarray hand_sides, + dpctl::tensor::usm_ndarray dependent_vals, const std::vector &depends); extern void init_gesv_dispatch_vector(void); diff --git a/dpnp/backend/extensions/lapack/lapack_py.cpp b/dpnp/backend/extensions/lapack/lapack_py.cpp index 5ae26d595cb8..feb96dc0772c 100644 --- a/dpnp/backend/extensions/lapack/lapack_py.cpp +++ b/dpnp/backend/extensions/lapack/lapack_py.cpp @@ -53,6 +53,7 @@ void init_dispatch_tables(void) PYBIND11_MODULE(_lapack_impl, m) { + // Register a custom LinAlgError exception in the dpnp.linalg submodule py::module_ linalg_module = py::module_::import("dpnp.linalg"); py::register_exception( linalg_module, "LinAlgError", PyExc_ValueError); @@ -76,9 +77,8 @@ PYBIND11_MODULE(_lapack_impl, m) m.def("_gesv", &lapack_ext::gesv, "Call `gesv` from OneMKL LAPACK library to return " - "solution to the system of linear equations with a square " - "coefficient matrix A" - "and multiple right-hand sides", - py::arg("sycl_queue"), py::arg("coeff_matrix"), py::arg("hand_sides"), - py::arg("depends") = py::list()); + "the solution of a system of linear equations with " + "a square coefficient matrix A and multiple dependent variables", + py::arg("sycl_queue"), py::arg("coeff_matrix"), + py::arg("dependent_vals"), py::arg("depends") = py::list()); } diff --git a/dpnp/backend/extensions/lapack/linalg_exceptions.hpp b/dpnp/backend/extensions/lapack/linalg_exceptions.hpp index 9963b2aa834e..083be22429c0 100644 --- a/dpnp/backend/extensions/lapack/linalg_exceptions.hpp +++ b/dpnp/backend/extensions/lapack/linalg_exceptions.hpp @@ -1,3 +1,28 @@ +//***************************************************************************** +// 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 diff --git a/dpnp/linalg/dpnp_iface_linalg.py b/dpnp/linalg/dpnp_iface_linalg.py index b9a0d889e046..3558010a8301 100644 --- a/dpnp/linalg/dpnp_iface_linalg.py +++ b/dpnp/linalg/dpnp_iface_linalg.py @@ -541,14 +541,14 @@ def solve(a, b): if a.ndim < 2: raise dpnp.linalg.LinAlgError( - f"{a.ndim}-dimensional array given. Array must be " - "at least two-dimensional" + f"{a.ndim}-dimensional array given. The input coefficient " + "array must be at least two-dimensional" ) m, n = a.shape[-2:] if m != n: raise dpnp.linalg.LinAlgError( - "Last 2 dimensions of the array must be square" + "Last 2 dimensions of the input coefficient array must be square" ) if not ( diff --git a/dpnp/linalg/dpnp_utils_linalg.py b/dpnp/linalg/dpnp_utils_linalg.py index d6bcd92fb788..a6c5ceadc8ab 100644 --- a/dpnp/linalg/dpnp_utils_linalg.py +++ b/dpnp/linalg/dpnp_utils_linalg.py @@ -173,7 +173,8 @@ def dpnp_solve(a, b): dpnp_solve(a, b) Return the the solution to the system of linear equations with - a square coefficient matrix `a` and multiple right-hand sides `b`. + a square coefficient matrix `a` and multiple dependent variables + array `b`. """ @@ -237,12 +238,13 @@ def dpnp_solve(a, b): for i in range(op_count): # oneMKL LAPACK assumes fortran-like array as input, so # allocate a memory with 'F' order for dpnp array of coefficient matrix - # and multiple right-hand sides + # and multiple dependent variables array coeff_vecs[i] = dpnp.empty_like(a[i], order="F", dtype=res_type) val_vecs[i] = dpnp.empty_like(b[i], order="F", dtype=res_type) - # use DPCTL tensor function to fill the array of coefficient matrix - # and multiple right-hand sides with content of input array + # use DPCTL tensor function to fill the coefficient matrix array + # and the array of multiple dependent variables with content + # from the input arrays a_ht_copy_ev[i], a_copy_ev = ti._copy_usm_ndarray_into_usm_ndarray( src=a_usm_arr[i], dst=coeff_vecs[i].get_array(), @@ -254,8 +256,9 @@ def dpnp_solve(a, b): sycl_queue=b.sycl_queue, ) - # call LAPACK extension function to get the solution of the system of linear - # equations with a portion of the coefficients square matrix + # Call the LAPACK extension function _gesv to solve the system of linear + # equations using a portion of the coefficient square matrix and a + # corresponding portion of the dependent variables array. ht_lapack_ev[i] = li._gesv( exec_q, coeff_vecs[i].get_array(), @@ -271,18 +274,20 @@ def dpnp_solve(a, b): # combine the list of solutions into a single array out_v = dpnp.array(val_vecs, order=b_order) if reshape: - # shape of the out_t must be equal to the shape of the right-hand sides + # shape of the out_v must be equal to the shape of the array of + # dependent variables out_v = out_v.reshape(orig_shape_b) return out_v else: # oneMKL LAPACK assumes fortran-like array as input, so # allocate a memory with 'F' order for dpnp array of coefficient matrix - # and multiple right-hand sides + # and multiple dependent variables a_f = dpnp.empty_like(a, order="F", dtype=res_type) b_f = dpnp.empty_like(b, order="F", dtype=res_type) - # use DPCTL tensor function to fill the array of coefficient matrix - # and multiple right-hand sides with content of input array + # use DPCTL tensor function to fill the coefficient matrix array + # and the array of multiple dependent variables with content + # from the input arrays a_ht_copy_ev, a_copy_ev = ti._copy_usm_ndarray_into_usm_ndarray( src=a_usm_arr, dst=a_f.get_array(), sycl_queue=a.sycl_queue ) @@ -290,14 +295,15 @@ def dpnp_solve(a, b): src=b_usm_arr, dst=b_f.get_array(), sycl_queue=b.sycl_queue ) - # call LAPACK extension function to get the solution of the system of linear - # equations with the coefficients square matrix + # Call the LAPACK extension function _gesv to solve the system of linear + # equations with the coefficient square matrix and the dependent variables array. lapack_ev = li._gesv( exec_q, a_f.get_array(), b_f.get_array(), [a_copy_ev, b_copy_ev] ) if b_order != "F": - # need to align order of the result of solutions with the right-hand sides + # need to align order of the result of solutions with the + # input array of multiple dependent variables out_v = dpnp.empty_like(b_f, order=b_order) ht_copy_out_ev, _ = ti._copy_usm_ndarray_into_usm_ndarray( src=b_f.get_array(), From 3f88dc5348f419f25d1a968c73e858b68e94017d Mon Sep 17 00:00:00 2001 From: Vladislav Perevezentsev Date: Mon, 16 Oct 2023 14:52:25 +0200 Subject: [PATCH 17/57] Simplify ThresholdType determination --- dpnp/backend/extensions/lapack/gesv.cpp | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/dpnp/backend/extensions/lapack/gesv.cpp b/dpnp/backend/extensions/lapack/gesv.cpp index 94b4f326eee1..52dd5e9fd2e2 100644 --- a/dpnp/backend/extensions/lapack/gesv.cpp +++ b/dpnp/backend/extensions/lapack/gesv.cpp @@ -114,13 +114,13 @@ static sycl::event gesv_impl(sycl::queue exec_q, exec_q.memcpy(&host_U, &a[(info - 1) * lda + info - 1], sizeof(T)) .wait(); - using ThresholdType = typename std::conditional< - std::is_same::value, float, - typename std::conditional< - std::is_same::value, double, - typename std::conditional< - std::is_same>::value, float, - double>::type>::type>::type; + using ThresholdType = std::conditional_t< + std::is_same_v || + std::is_same_v>, + float, + std::conditional_t || + std::is_same_v>, + double, void>>; const auto threshold = std::numeric_limits::epsilon() * 100; From 22e87343420744b2963739cb13626c084f9c34ea Mon Sep 17 00:00:00 2001 From: Vladislav Perevezentsev Date: Tue, 17 Oct 2023 10:58:51 +0200 Subject: [PATCH 18/57] Small changes to the docstrings --- dpnp/backend/extensions/lapack/types_matrix.hpp | 5 +++-- dpnp/linalg/dpnp_utils_linalg.py | 2 +- 2 files changed, 4 insertions(+), 3 deletions(-) diff --git a/dpnp/backend/extensions/lapack/types_matrix.hpp b/dpnp/backend/extensions/lapack/types_matrix.hpp index 9d0a55c41ccf..124774982e8d 100644 --- a/dpnp/backend/extensions/lapack/types_matrix.hpp +++ b/dpnp/backend/extensions/lapack/types_matrix.hpp @@ -86,8 +86,9 @@ struct SyevdTypePairSupportFactory * MKL LAPACK library provides support in oneapi::mkl::lapack::gesv * function. * - * @tparam T Type of array containing input matrix A and an output arrays with - * coefficient matrix and hand sides. + * @tparam T Type of array containing the coefficient matrix A and + * the array of multiple dependent variables. Upon execution, the array of + * multiple dependent variables will be overwritten with the solution. */ template struct GesvTypePairSupportFactory diff --git a/dpnp/linalg/dpnp_utils_linalg.py b/dpnp/linalg/dpnp_utils_linalg.py index efb481b30a3e..6563ed7df544 100644 --- a/dpnp/linalg/dpnp_utils_linalg.py +++ b/dpnp/linalg/dpnp_utils_linalg.py @@ -172,7 +172,7 @@ def dpnp_solve(a, b): """ dpnp_solve(a, b) - Return the the solution to the system of linear equations with + Return the solution to the system of linear equations with a square coefficient matrix `a` and multiple dependent variables array `b`. From 87cda0be528cf285b6b142f16f52c5386348cc72 Mon Sep 17 00:00:00 2001 From: Vladislav Perevezentsev Date: Tue, 17 Oct 2023 17:16:02 +0200 Subject: [PATCH 19/57] Remove if op_count due to unreachable --- dpnp/linalg/dpnp_utils_linalg.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/dpnp/linalg/dpnp_utils_linalg.py b/dpnp/linalg/dpnp_utils_linalg.py index 6563ed7df544..e360c0ab22eb 100644 --- a/dpnp/linalg/dpnp_utils_linalg.py +++ b/dpnp/linalg/dpnp_utils_linalg.py @@ -226,8 +226,6 @@ def dpnp_solve(a, b): reshape = True op_count = a.shape[0] - if op_count == 0: - return dpnp.empty_like(b, dtype=res_type) coeff_vecs = [None] * op_count val_vecs = [None] * op_count From 76e035d4c2f8fd27f2c5656979587389295289e3 Mon Sep 17 00:00:00 2001 From: Vladislav Perevezentsev Date: Tue, 17 Oct 2023 17:27:10 +0200 Subject: [PATCH 20/57] Improve test coverage --- tests/test_linalg.py | 26 ++++++++++++++++++++++++++ 1 file changed, 26 insertions(+) diff --git a/tests/test_linalg.py b/tests/test_linalg.py index 8392526817b3..a223f449cf38 100644 --- a/tests/test_linalg.py +++ b/tests/test_linalg.py @@ -545,3 +545,29 @@ def test_solve_strides(self): expected = numpy.linalg.solve(a_np[::-2, ::-2], b_np[::-2]) result = inp.linalg.solve(a_dp[::-2, ::-2], b_dp[::-2]) assert_allclose(expected, result, rtol=1e-05) + + def test_solve_errors(self): + # different type + a_dp = inp.array([[1, 0.5], [0.5, 1]], dtype="float32") + b_dp = inp.array(a_dp, dtype="float32") + a_dp_int = a_dp.astype("int32") + assert_raises(ValueError, inp.linalg.solve, a_dp_int, b_dp) + + # diffetent queue + a_queue = dpctl.SyclQueue() + b_queue = dpctl.SyclQueue() + a_dp_q = inp.array(a_dp, sycl_queue=a_queue) + b_dp_q = inp.array(b_dp, sycl_queue=b_queue) + assert_raises(ValueError, inp.linalg.solve, a_dp_q, b_dp_q) + + # unsupported type + a_np = inp.asnumpy(a_dp) + b_np = inp.asnumpy(b_dp) + assert_raises(TypeError, inp.linalg.solve, a_np, b_dp) + assert_raises(TypeError, inp.linalg.solve, a_dp, b_np) + + # a.ndim < 2 + a_dp_ndim_1 = a_dp.flatten() + assert_raises( + inp.linalg.LinAlgError, inp.linalg.solve, a_dp_ndim_1, b_dp + ) From 1bfc81a89c61c2f5e0ac8c46741f3802ca555083 Mon Sep 17 00:00:00 2001 From: Vladislav Perevezentsev Date: Wed, 18 Oct 2023 00:02:50 +0200 Subject: [PATCH 21/57] Impl dtype dispatching with linalg_common_type for dpnp.linalg.solve --- dpnp/linalg/dpnp_utils_linalg.py | 75 ++++++++++++++++++++++---------- 1 file changed, 52 insertions(+), 23 deletions(-) diff --git a/dpnp/linalg/dpnp_utils_linalg.py b/dpnp/linalg/dpnp_utils_linalg.py index e360c0ab22eb..1aef92399e70 100644 --- a/dpnp/linalg/dpnp_utils_linalg.py +++ b/dpnp/linalg/dpnp_utils_linalg.py @@ -34,12 +34,53 @@ import dpnp import dpnp.backend.extensions.lapack._lapack_impl as li -__all__ = ["dpnp_eigh"] +__all__ = ["dpnp_eigh", "dpnp_solve"] _jobz = {"N": 0, "V": 1} _upper_lower = {"U": 0, "L": 1} +def linalg_common_type(*arrays): + """ + linalg_common_type(*arrays) + + Common type for linalg + + This function determines the common data type for linalg operations. + It's designed to be similar in logic to `numpy.linalg.linalg._commonType`. + + Key differences from `numpy.common_type`: + - It accepts ``bool_`` arrays. + - It doesn't consider ``float16`` arrays since they are not supported. + - The default floating-point data type is determined by the capabilities of the device, + as indicated by `dpnp.default_float_type()`. + + Args: + *arrays (dpnp.ndarray): Input arrays. + + Returns: + dtype_common (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). + + """ + + dtypes = [arr.dtype for arr in arrays] + + default = dpnp.default_float_type().name + dtype_common = _common_type_internal(default, *dtypes) + + return dtype_common, dtype_common + + +def _common_type_internal(default_dtype, *dtypes): + inexact_dtypes = [ + dtype if dtype.kind in "fc" else default_dtype for dtype in dtypes + ] + return dpnp.result_type(*inexact_dtypes) + + def dpnp_eigh(a, UPLO): """ dpnp_eigh(a, UPLO) @@ -182,9 +223,8 @@ def dpnp_solve(a, b): b_usm_arr = dpnp.get_usm_ndarray(b) b_order = "C" if b.flags.c_contiguous else "F" - - if a.dtype != b.dtype: - raise ValueError("a and b must be of the same type") + a_shape = a.shape + b_shape = b.shape exec_q = dpctl.utils.get_execution_queue((a.sycl_queue, b.sycl_queue)) if exec_q is None: @@ -193,33 +233,22 @@ def dpnp_solve(a, b): "from input arguments." ) - if dpnp.issubdtype(a.dtype, dpnp.floating): - res_type = ( - a.dtype if exec_q.sycl_device.has_aspect_fp64 else dpnp.float32 - ) - elif dpnp.issubdtype(a.dtype, dpnp.complexfloating): - res_type = ( - a.dtype if exec_q.sycl_device.has_aspect_fp64 else dpnp.complex64 - ) - else: - res_type = ( - dpnp.float64 if exec_q.sycl_device.has_aspect_fp64 else dpnp.float32 - ) + dtype, res_type = linalg_common_type(a, b) if b.size == 0: - return dpnp.empty(b.shape, dtype=res_type) + return dpnp.empty(b_shape, dtype=res_type) if a.ndim > 2: reshape = False - orig_shape_b = b.shape + orig_shape_b = b_shape if a.ndim > 3: # get 3d input arrays by reshape if a.ndim == b.ndim: - b = b.reshape(prod(b.shape[:-2]), b.shape[-2], b.shape[-1]) + b = b.reshape(prod(b_shape[:-2]), b_shape[-2], b_shape[-1]) else: - b = b.reshape(prod(b.shape[:-1]), b.shape[-1]) + b = b.reshape(prod(b_shape[:-1]), b_shape[-1]) - a = a.reshape(prod(a.shape[:-2]), a.shape[-2], a.shape[-1]) + a = a.reshape(prod(a_shape[:-2]), a_shape[-2], a_shape[-1]) a_usm_arr = dpnp.get_usm_ndarray(a) b_usm_arr = dpnp.get_usm_ndarray(b) @@ -237,7 +266,7 @@ def dpnp_solve(a, b): # oneMKL LAPACK assumes fortran-like array as input, so # allocate a memory with 'F' order for dpnp array of coefficient matrix # and multiple dependent variables array - coeff_vecs[i] = dpnp.empty_like(a[i], order="F", dtype=res_type) + coeff_vecs[i] = dpnp.empty_like(a[i], order="F", dtype=dtype) val_vecs[i] = dpnp.empty_like(b[i], order="F", dtype=res_type) # use DPCTL tensor function to fill the coefficient matrix array @@ -280,7 +309,7 @@ def dpnp_solve(a, b): # oneMKL LAPACK assumes fortran-like array as input, so # allocate a memory with 'F' order for dpnp array of coefficient matrix # and multiple dependent variables - a_f = dpnp.empty_like(a, order="F", dtype=res_type) + a_f = dpnp.empty_like(a, order="F", dtype=dtype) b_f = dpnp.empty_like(b, order="F", dtype=res_type) # use DPCTL tensor function to fill the coefficient matrix array From c7b284b81bea086f92804166715a197fede443ef Mon Sep 17 00:00:00 2001 From: Vladislav Perevezentsev Date: Wed, 18 Oct 2023 00:03:40 +0200 Subject: [PATCH 22/57] Add a new test_solve_diff_type --- tests/test_linalg.py | 28 ++++++++++++++++++++++------ 1 file changed, 22 insertions(+), 6 deletions(-) diff --git a/tests/test_linalg.py b/tests/test_linalg.py index a223f449cf38..396e60ae9883 100644 --- a/tests/test_linalg.py +++ b/tests/test_linalg.py @@ -5,7 +5,7 @@ import dpnp as inp -from .helper import get_all_dtypes, get_complex_dtypes, has_support_aspect64 +from .helper import get_all_dtypes, has_support_aspect64 def vvsort(val, vec, size, xp): @@ -513,14 +513,33 @@ def test_svd(type, shape): class TestSolve: @pytest.mark.parametrize("dtype", get_all_dtypes(no_bool=True)) def test_solve(self, dtype): - a_np = numpy.array([[1, 0.5], [0.5, 1]]) - a_dp = inp.array(a_np, dtype=dtype) + a_np = numpy.array([[1, 0.5], [0.5, 1]], dtype=dtype) + a_dp = inp.array(a_np) expected = numpy.linalg.solve(a_np, a_np) result = inp.linalg.solve(a_dp, a_dp) assert_allclose(expected, result, rtol=1e-06) + @pytest.mark.parametrize("a_dtype", get_all_dtypes(no_bool=True)) + @pytest.mark.parametrize("b_dtype", get_all_dtypes(no_bool=True)) + def test_solve_diff_type(self, a_dtype, b_dtype): + a_np = numpy.array([[1, 2], [3, -5]], dtype=a_dtype) + b_np = numpy.array([4, 1], dtype=b_dtype) + + a_dp = inp.array(a_np) + b_dp = inp.array(b_np) + + expected = numpy.linalg.solve(a_np, b_np) + result = inp.linalg.solve(a_dp, b_dp) + + assert_allclose(expected, result, rtol=1e-06) + + if has_support_aspect64(): + assert expected.dtype == result.dtype + else: + assert expected.dtype.kind == result.dtype.kind + def test_solve_strides(self): a_np = numpy.array( [ @@ -547,11 +566,8 @@ def test_solve_strides(self): assert_allclose(expected, result, rtol=1e-05) def test_solve_errors(self): - # different type a_dp = inp.array([[1, 0.5], [0.5, 1]], dtype="float32") b_dp = inp.array(a_dp, dtype="float32") - a_dp_int = a_dp.astype("int32") - assert_raises(ValueError, inp.linalg.solve, a_dp_int, b_dp) # diffetent queue a_queue = dpctl.SyclQueue() From e99d37c6e6db12411f082640f15915e5d36d71bc Mon Sep 17 00:00:00 2001 From: Vladislav Perevezentsev Date: Thu, 26 Oct 2023 11:48:03 +0200 Subject: [PATCH 23/57] Add a common_helpers.hpp file --- .../extensions/lapack/common_helpers.hpp | 55 +++++++++++++++++++ dpnp/backend/extensions/lapack/gesv.cpp | 9 +-- 2 files changed, 57 insertions(+), 7 deletions(-) create mode 100644 dpnp/backend/extensions/lapack/common_helpers.hpp diff --git a/dpnp/backend/extensions/lapack/common_helpers.hpp b/dpnp/backend/extensions/lapack/common_helpers.hpp new file mode 100644 index 000000000000..2f3815320cad --- /dev/null +++ b/dpnp/backend/extensions/lapack/common_helpers.hpp @@ -0,0 +1,55 @@ +//***************************************************************************** +// 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 + +namespace dpnp +{ +namespace backend +{ +namespace ext +{ +namespace lapack +{ +namespace helper +{ +template +struct value_type_of +{ + using type = T; +}; + +template +struct value_type_of> +{ + using type = T; +}; +} // namespace helper +} // namespace lapack +} // namespace ext +} // namespace backend +} // namespace dpnp diff --git a/dpnp/backend/extensions/lapack/gesv.cpp b/dpnp/backend/extensions/lapack/gesv.cpp index 52dd5e9fd2e2..763b602411a4 100644 --- a/dpnp/backend/extensions/lapack/gesv.cpp +++ b/dpnp/backend/extensions/lapack/gesv.cpp @@ -29,6 +29,7 @@ #include "utils/memory_overlap.hpp" #include "utils/type_utils.hpp" +#include "common_helpers.hpp" #include "gesv.hpp" #include "linalg_exceptions.hpp" #include "types_matrix.hpp" @@ -114,13 +115,7 @@ static sycl::event gesv_impl(sycl::queue exec_q, exec_q.memcpy(&host_U, &a[(info - 1) * lda + info - 1], sizeof(T)) .wait(); - using ThresholdType = std::conditional_t< - std::is_same_v || - std::is_same_v>, - float, - std::conditional_t || - std::is_same_v>, - double, void>>; + using ThresholdType = typename helper::value_type_of::type; const auto threshold = std::numeric_limits::epsilon() * 100; From ec1e966e994d52b78bbe33e39882c37b559f40a1 Mon Sep 17 00:00:00 2001 From: Vladislav Perevezentsev Date: Thu, 26 Oct 2023 11:55:25 +0200 Subject: [PATCH 24/57] Use bool flag for sycl exception --- dpnp/backend/extensions/lapack/gesv.cpp | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/dpnp/backend/extensions/lapack/gesv.cpp b/dpnp/backend/extensions/lapack/gesv.cpp index 763b602411a4..7165e1aa2de5 100644 --- a/dpnp/backend/extensions/lapack/gesv.cpp +++ b/dpnp/backend/extensions/lapack/gesv.cpp @@ -84,6 +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; sycl::event gesv_event; try { @@ -142,10 +143,10 @@ static sycl::event gesv_impl(sycl::queue exec_q, } catch (sycl::exception const &e) { error_msg << "Unexpected SYCL exception caught during gesv() call:\n" << e.what(); - info = -11; + sycl_exception_caught = true; } - if (info != 0) // an unexected error occurs + if (info != 0 || sycl_exception_caught) // an unexected error occurs { if (scratchpad != nullptr) { sycl::free(scratchpad, exec_q); From 6885aa3f5fb953a2453deaa9eaec23b39282e49a Mon Sep 17 00:00:00 2001 From: Vladislav Perevezentsev Date: Thu, 26 Oct 2023 12:14:16 +0200 Subject: [PATCH 25/57] Refactor memory management for ipiv in gesv_impl --- dpnp/backend/extensions/lapack/gesv.cpp | 20 ++++++++++++-------- 1 file changed, 12 insertions(+), 8 deletions(-) diff --git a/dpnp/backend/extensions/lapack/gesv.cpp b/dpnp/backend/extensions/lapack/gesv.cpp index 7165e1aa2de5..56e54f593215 100644 --- a/dpnp/backend/extensions/lapack/gesv.cpp +++ b/dpnp/backend/extensions/lapack/gesv.cpp @@ -53,7 +53,6 @@ typedef sycl::event (*gesv_impl_fn_ptr_t)(sycl::queue, const std::int64_t, char *, std::int64_t, - std::int64_t *, char *, std::int64_t, std::vector &, @@ -67,7 +66,6 @@ static sycl::event gesv_impl(sycl::queue exec_q, const std::int64_t nrhs, char *in_a, std::int64_t lda, - std::int64_t *ipiv, char *in_b, std::int64_t ldb, std::vector &host_task_events, @@ -82,6 +80,8 @@ static sycl::event gesv_impl(sycl::queue exec_q, mkl_lapack::gesv_scratchpad_size(exec_q, n, nrhs, lda, ldb); T *scratchpad = nullptr; + std::int64_t *ipiv = nullptr; + std::stringstream error_msg; std::int64_t info = 0; bool sycl_exception_caught = false; @@ -89,6 +89,7 @@ static sycl::event gesv_impl(sycl::queue exec_q, sycl::event gesv_event; try { scratchpad = sycl::malloc_device(scratchpad_size, exec_q); + ipiv = sycl::malloc_device(n, exec_q); gesv_event = mkl_lapack::gesv( exec_q, @@ -151,13 +152,19 @@ static sycl::event gesv_impl(sycl::queue exec_q, if (scratchpad != nullptr) { sycl::free(scratchpad, exec_q); } + if (ipiv != nullptr) { + sycl::free(ipiv, exec_q); + } throw std::runtime_error(error_msg.str()); } sycl::event clean_up_event = exec_q.submit([&](sycl::handler &cgh) { cgh.depends_on(gesv_event); auto ctx = exec_q.get_context(); - cgh.host_task([ctx, scratchpad]() { sycl::free(scratchpad, ctx); }); + cgh.host_task([ctx, scratchpad, ipiv]() { + sycl::free(scratchpad, ctx); + sycl::free(ipiv, ctx); + }); }); host_task_events.push_back(clean_up_event); @@ -238,13 +245,10 @@ sycl::event gesv(sycl::queue exec_q, const std::int64_t lda = std::max(1UL, n); const std::int64_t ldb = std::max(1UL, m); - std::vector ipiv(n); - std::int64_t *d_ipiv = sycl::malloc_device(n, exec_q); - std::vector host_task_events; sycl::event gesv_ev = - gesv_fn(exec_q, n, nrhs, coeff_matrix_data, lda, d_ipiv, - dependent_vals_data, ldb, host_task_events, depends); + gesv_fn(exec_q, n, nrhs, coeff_matrix_data, lda, dependent_vals_data, + ldb, host_task_events, depends); return gesv_ev; } From 56c920f8ccac486b42f918eb5c33504960a95ca9 Mon Sep 17 00:00:00 2001 From: Vladislav Perevezentsev Date: Fri, 27 Oct 2023 13:19:04 +0200 Subject: [PATCH 26/57] Rename linalg_common_type to _common_type and change the number of types it returns --- dpnp/linalg/dpnp_utils_linalg.py | 13 ++++++------- 1 file changed, 6 insertions(+), 7 deletions(-) diff --git a/dpnp/linalg/dpnp_utils_linalg.py b/dpnp/linalg/dpnp_utils_linalg.py index 1aef92399e70..81d55937329b 100644 --- a/dpnp/linalg/dpnp_utils_linalg.py +++ b/dpnp/linalg/dpnp_utils_linalg.py @@ -40,9 +40,9 @@ _upper_lower = {"U": 0, "L": 1} -def linalg_common_type(*arrays): +def _common_type(*arrays): """ - linalg_common_type(*arrays) + _common_type(*arrays) Common type for linalg @@ -71,7 +71,7 @@ def linalg_common_type(*arrays): default = dpnp.default_float_type().name dtype_common = _common_type_internal(default, *dtypes) - return dtype_common, dtype_common + return dtype_common def _common_type_internal(default_dtype, *dtypes): @@ -233,8 +233,7 @@ def dpnp_solve(a, b): "from input arguments." ) - dtype, res_type = linalg_common_type(a, b) - + res_type = _common_type(a, b) if b.size == 0: return dpnp.empty(b_shape, dtype=res_type) @@ -266,7 +265,7 @@ def dpnp_solve(a, b): # oneMKL LAPACK assumes fortran-like array as input, so # allocate a memory with 'F' order for dpnp array of coefficient matrix # and multiple dependent variables array - coeff_vecs[i] = dpnp.empty_like(a[i], order="F", dtype=dtype) + coeff_vecs[i] = dpnp.empty_like(a[i], order="F", dtype=res_type) val_vecs[i] = dpnp.empty_like(b[i], order="F", dtype=res_type) # use DPCTL tensor function to fill the coefficient matrix array @@ -309,7 +308,7 @@ def dpnp_solve(a, b): # oneMKL LAPACK assumes fortran-like array as input, so # allocate a memory with 'F' order for dpnp array of coefficient matrix # and multiple dependent variables - a_f = dpnp.empty_like(a, order="F", dtype=dtype) + a_f = dpnp.empty_like(a, order="F", dtype=res_type) b_f = dpnp.empty_like(b, order="F", dtype=res_type) # use DPCTL tensor function to fill the coefficient matrix array From 07b1418cbbe7be9bbb689635b208940e88212a4c Mon Sep 17 00:00:00 2001 From: Vladislav Perevezentsev Date: Fri, 27 Oct 2023 13:20:49 +0200 Subject: [PATCH 27/57] Address the remarks --- dpnp/backend/extensions/lapack/lapack_py.cpp | 14 ++--- .../extensions/lapack/types_matrix.hpp | 53 +++++++++---------- dpnp/linalg/dpnp_iface_linalg.py | 1 + 3 files changed, 34 insertions(+), 34 deletions(-) diff --git a/dpnp/backend/extensions/lapack/lapack_py.cpp b/dpnp/backend/extensions/lapack/lapack_py.cpp index feb96dc0772c..4273c90697a4 100644 --- a/dpnp/backend/extensions/lapack/lapack_py.cpp +++ b/dpnp/backend/extensions/lapack/lapack_py.cpp @@ -61,6 +61,13 @@ PYBIND11_MODULE(_lapack_impl, m) init_dispatch_vectors(); init_dispatch_tables(); + m.def("_gesv", &lapack_ext::gesv, + "Call `gesv` from OneMKL LAPACK library to return " + "the solution of a system of linear equations with " + "a square coefficient matrix A and multiple dependent variables", + py::arg("sycl_queue"), py::arg("coeff_matrix"), + py::arg("dependent_vals"), 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", @@ -74,11 +81,4 @@ PYBIND11_MODULE(_lapack_impl, m) py::arg("sycl_queue"), py::arg("jobz"), py::arg("upper_lower"), py::arg("eig_vecs"), py::arg("eig_vals"), py::arg("depends") = py::list()); - - m.def("_gesv", &lapack_ext::gesv, - "Call `gesv` from OneMKL LAPACK library to return " - "the solution of a system of linear equations with " - "a square coefficient matrix A and multiple dependent variables", - py::arg("sycl_queue"), py::arg("coeff_matrix"), - py::arg("dependent_vals"), py::arg("depends") = py::list()); } diff --git a/dpnp/backend/extensions/lapack/types_matrix.hpp b/dpnp/backend/extensions/lapack/types_matrix.hpp index 124774982e8d..60521cb75a3c 100644 --- a/dpnp/backend/extensions/lapack/types_matrix.hpp +++ b/dpnp/backend/extensions/lapack/types_matrix.hpp @@ -43,6 +43,32 @@ namespace lapack { namespace types { +/** + * @brief A factory to define pairs of supported types for which + * MKL LAPACK library provides support in oneapi::mkl::lapack::gesv + * function. + * + * @tparam T Type of array containing the coefficient matrix A and + * the array of multiple dependent variables. Upon execution, the array of + * multiple dependent variables will be overwritten with the solution. + */ +template +struct GesvTypePairSupportFactory +{ + 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 @@ -80,33 +106,6 @@ struct SyevdTypePairSupportFactory // 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::gesv - * function. - * - * @tparam T Type of array containing the coefficient matrix A and - * the array of multiple dependent variables. Upon execution, the array of - * multiple dependent variables will be overwritten with the solution. - */ -template -struct GesvTypePairSupportFactory -{ - 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; -}; } // namespace types } // namespace lapack } // namespace ext diff --git a/dpnp/linalg/dpnp_iface_linalg.py b/dpnp/linalg/dpnp_iface_linalg.py index 25383eb9a4e5..5162a4545303 100644 --- a/dpnp/linalg/dpnp_iface_linalg.py +++ b/dpnp/linalg/dpnp_iface_linalg.py @@ -533,6 +533,7 @@ def solve(a, b): array([-1., 1.]) """ + if not dpnp.is_supported_array_type(a): raise TypeError( "An array must be any of supported type, but got {}".format(type(a)) From 0336c002300977df382dd4c7d96a96ecec1a9806 Mon Sep 17 00:00:00 2001 From: Vladislav Perevezentsev Date: Thu, 2 Nov 2023 11:02:41 +0100 Subject: [PATCH 28/57] Remove the use of prod to get 3d array and rename op_count to batch_size --- dpnp/linalg/dpnp_utils_linalg.py | 25 ++++++++++++------------- 1 file changed, 12 insertions(+), 13 deletions(-) diff --git a/dpnp/linalg/dpnp_utils_linalg.py b/dpnp/linalg/dpnp_utils_linalg.py index 81d55937329b..87f40f3d6ab7 100644 --- a/dpnp/linalg/dpnp_utils_linalg.py +++ b/dpnp/linalg/dpnp_utils_linalg.py @@ -29,7 +29,6 @@ import dpctl import dpctl.tensor._tensor_impl as ti -from numpy import prod import dpnp import dpnp.backend.extensions.lapack._lapack_impl as li @@ -243,25 +242,25 @@ def dpnp_solve(a, b): if a.ndim > 3: # get 3d input arrays by reshape if a.ndim == b.ndim: - b = b.reshape(prod(b_shape[:-2]), b_shape[-2], b_shape[-1]) + b = b.reshape(-1, b_shape[-2], b_shape[-1]) else: - b = b.reshape(prod(b_shape[:-1]), b_shape[-1]) + b = b.reshape(-1, b_shape[-1]) - a = a.reshape(prod(a_shape[:-2]), a_shape[-2], a_shape[-1]) + a = a.reshape(-1, a_shape[-2], a_shape[-1]) a_usm_arr = dpnp.get_usm_ndarray(a) b_usm_arr = dpnp.get_usm_ndarray(b) reshape = True - op_count = a.shape[0] + batch_size = a.shape[0] - coeff_vecs = [None] * op_count - val_vecs = [None] * op_count - a_ht_copy_ev = [None] * op_count - b_ht_copy_ev = [None] * op_count - ht_lapack_ev = [None] * op_count + coeff_vecs = [None] * batch_size + val_vecs = [None] * batch_size + a_ht_copy_ev = [None] * batch_size + b_ht_copy_ev = [None] * batch_size + ht_lapack_ev = [None] * batch_size - for i in range(op_count): + for i in range(batch_size): # oneMKL LAPACK assumes fortran-like array as input, so # allocate a memory with 'F' order for dpnp array of coefficient matrix # and multiple dependent variables array @@ -292,7 +291,7 @@ def dpnp_solve(a, b): depends=[a_copy_ev, b_copy_ev], ) - for i in range(op_count): + for i in range(batch_size): ht_lapack_ev[i].wait() b_ht_copy_ev[i].wait() a_ht_copy_ev[i].wait() @@ -334,7 +333,7 @@ def dpnp_solve(a, b): ht_copy_out_ev, _ = ti._copy_usm_ndarray_into_usm_ndarray( src=b_f.get_array(), dst=out_v.get_array(), - sycl_queue=b.sycl_queue, + sycl_queue=exec_q, depends=[lapack_ev], ) ht_copy_out_ev.wait() From 31f6f10ca1031af84eda1a58ad58d12b151bad50 Mon Sep 17 00:00:00 2001 From: Vladislav Perevezentsev Date: Thu, 2 Nov 2023 11:23:54 +0100 Subject: [PATCH 29/57] gesv returns pair of events and uses dpctl.utils.keep_args_alive --- dpnp/backend/extensions/lapack/gesv.cpp | 14 +++++++++----- dpnp/backend/extensions/lapack/gesv.hpp | 9 +++++---- dpnp/linalg/dpnp_utils_linalg.py | 6 +++--- 3 files changed, 17 insertions(+), 12 deletions(-) diff --git a/dpnp/backend/extensions/lapack/gesv.cpp b/dpnp/backend/extensions/lapack/gesv.cpp index 56e54f593215..7ac518f3ed26 100644 --- a/dpnp/backend/extensions/lapack/gesv.cpp +++ b/dpnp/backend/extensions/lapack/gesv.cpp @@ -171,10 +171,11 @@ static sycl::event gesv_impl(sycl::queue exec_q, return gesv_event; } -sycl::event gesv(sycl::queue exec_q, - dpctl::tensor::usm_ndarray coeff_matrix, - dpctl::tensor::usm_ndarray dependent_vals, - const std::vector &depends) +std::pair + gesv(sycl::queue exec_q, + dpctl::tensor::usm_ndarray coeff_matrix, + dpctl::tensor::usm_ndarray dependent_vals, + const std::vector &depends) { const int coeff_matrix_nd = coeff_matrix.get_ndim(); const int dependent_vals_nd = dependent_vals.get_ndim(); @@ -250,7 +251,10 @@ sycl::event gesv(sycl::queue exec_q, gesv_fn(exec_q, n, nrhs, coeff_matrix_data, lda, dependent_vals_data, ldb, host_task_events, depends); - return gesv_ev; + sycl::event args_ev = dpctl::utils::keep_args_alive( + exec_q, {coeff_matrix, dependent_vals}, host_task_events); + + return std::make_pair(args_ev, gesv_ev); } template diff --git a/dpnp/backend/extensions/lapack/gesv.hpp b/dpnp/backend/extensions/lapack/gesv.hpp index 3b9a099e9795..24ac0d2e5be1 100644 --- a/dpnp/backend/extensions/lapack/gesv.hpp +++ b/dpnp/backend/extensions/lapack/gesv.hpp @@ -38,10 +38,11 @@ namespace ext { namespace lapack { -extern sycl::event gesv(sycl::queue exec_q, - dpctl::tensor::usm_ndarray coeff_matrix, - dpctl::tensor::usm_ndarray dependent_vals, - const std::vector &depends); +extern std::pair + gesv(sycl::queue exec_q, + dpctl::tensor::usm_ndarray coeff_matrix, + dpctl::tensor::usm_ndarray dependent_vals, + const std::vector &depends); extern void init_gesv_dispatch_vector(void); } // namespace lapack diff --git a/dpnp/linalg/dpnp_utils_linalg.py b/dpnp/linalg/dpnp_utils_linalg.py index 87f40f3d6ab7..c3f47f1e7939 100644 --- a/dpnp/linalg/dpnp_utils_linalg.py +++ b/dpnp/linalg/dpnp_utils_linalg.py @@ -284,7 +284,7 @@ def dpnp_solve(a, b): # Call the LAPACK extension function _gesv to solve the system of linear # equations using a portion of the coefficient square matrix and a # corresponding portion of the dependent variables array. - ht_lapack_ev[i] = li._gesv( + ht_lapack_ev[i], _ = li._gesv( exec_q, coeff_vecs[i].get_array(), val_vecs[i].get_array(), @@ -322,7 +322,7 @@ def dpnp_solve(a, b): # Call the LAPACK extension function _gesv to solve the system of linear # equations with the coefficient square matrix and the dependent variables array. - lapack_ev = li._gesv( + ht_lapack_ev, lapack_ev = li._gesv( exec_q, a_f.get_array(), b_f.get_array(), [a_copy_ev, b_copy_ev] ) @@ -340,7 +340,7 @@ def dpnp_solve(a, b): else: out_v = b_f - lapack_ev.wait() + ht_lapack_ev.wait() b_ht_copy_ev.wait() a_ht_copy_ev.wait() From d3717a68ded7acd7bf1784fbec7d045480c02abd Mon Sep 17 00:00:00 2001 From: Vladislav Perevezentsev Date: Thu, 2 Nov 2023 13:13:21 +0100 Subject: [PATCH 30/57] Return the use prod to get 3d arrays --- dpnp/linalg/dpnp_utils_linalg.py | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/dpnp/linalg/dpnp_utils_linalg.py b/dpnp/linalg/dpnp_utils_linalg.py index c3f47f1e7939..7739223b8646 100644 --- a/dpnp/linalg/dpnp_utils_linalg.py +++ b/dpnp/linalg/dpnp_utils_linalg.py @@ -29,6 +29,7 @@ import dpctl import dpctl.tensor._tensor_impl as ti +from numpy import prod import dpnp import dpnp.backend.extensions.lapack._lapack_impl as li @@ -242,11 +243,11 @@ def dpnp_solve(a, b): if a.ndim > 3: # get 3d input arrays by reshape if a.ndim == b.ndim: - b = b.reshape(-1, b_shape[-2], b_shape[-1]) + b = b.reshape(prod(b_shape[:-2]), b_shape[-2], b_shape[-1]) else: - b = b.reshape(-1, b_shape[-1]) + b = b.reshape(prod(b_shape[:-1]), b_shape[-1]) - a = a.reshape(-1, a_shape[-2], a_shape[-1]) + a = a.reshape(prod(a_shape[:-2]), a_shape[-2], a_shape[-1]) a_usm_arr = dpnp.get_usm_ndarray(a) b_usm_arr = dpnp.get_usm_ndarray(b) From 515df4e3f8b2f5a0de992a4de9479623ba89ce68 Mon Sep 17 00:00:00 2001 From: Vladislav Perevezentsev Date: Mon, 20 Nov 2023 16:10:16 +0100 Subject: [PATCH 31/57] Add test_solve_singular_matrix in TestSolve --- tests/test_linalg.py | 29 +++++++++++++++++++++++++++++ 1 file changed, 29 insertions(+) diff --git a/tests/test_linalg.py b/tests/test_linalg.py index 396e60ae9883..8162bd4081a8 100644 --- a/tests/test_linalg.py +++ b/tests/test_linalg.py @@ -565,6 +565,35 @@ def test_solve_strides(self): result = inp.linalg.solve(a_dp[::-2, ::-2], b_dp[::-2]) assert_allclose(expected, result, rtol=1e-05) + @pytest.mark.parametrize( + "matrix, vector", + [ + ([[1, 2], [2, 4]], [1, 2]), + ([[0, 0], [0, 0]], [0, 0]), + ([[1, 1], [1, 1]], [2, 2]), + ([[2, 4], [1, 2]], [3, 1.5]), + ([[1, 2], [0, 0]], [3, 0]), + ([[1, 0], [2, 0]], [3, 4]), + ], + ids=[ + "Linearly dependent rows", + "Zero matrix", + "Identical rows", + "Linearly dependent columns", + "Zero row", + "Zero column", + ], + ) + def test_solve_singular_matrix(self, matrix, vector): + a_np = numpy.array(matrix, dtype="float32") + b_np = numpy.array(vector, dtype="float32") + + a_dp = inp.array(a_np) + b_dp = inp.array(b_np) + + assert_raises(numpy.linalg.LinAlgError, numpy.linalg.solve, a_np, b_np) + assert_raises(inp.linalg.LinAlgError, inp.linalg.solve, a_dp, b_dp) + def test_solve_errors(self): a_dp = inp.array([[1, 0.5], [0.5, 1]], dtype="float32") b_dp = inp.array(a_dp, dtype="float32") From f992b0beb13a6319e464075d893486cdeea8bfc9 Mon Sep 17 00:00:00 2001 From: Vladislav Perevezentsev Date: Mon, 20 Nov 2023 16:26:09 +0100 Subject: [PATCH 32/57] Adress the remarks --- dpnp/backend/extensions/lapack/gesv.cpp | 2 +- dpnp/backend/extensions/lapack/lapack_py.cpp | 2 +- dpnp/linalg/dpnp_utils_linalg.py | 3 +-- tests/test_linalg.py | 8 ++------ 4 files changed, 5 insertions(+), 10 deletions(-) diff --git a/dpnp/backend/extensions/lapack/gesv.cpp b/dpnp/backend/extensions/lapack/gesv.cpp index 7ac518f3ed26..f536564fb92f 100644 --- a/dpnp/backend/extensions/lapack/gesv.cpp +++ b/dpnp/backend/extensions/lapack/gesv.cpp @@ -147,7 +147,7 @@ static sycl::event gesv_impl(sycl::queue exec_q, sycl_exception_caught = true; } - if (info != 0 || sycl_exception_caught) // an unexected error occurs + if (info != 0 || sycl_exception_caught) // an unexpected error occurs { if (scratchpad != nullptr) { sycl::free(scratchpad, exec_q); diff --git a/dpnp/backend/extensions/lapack/lapack_py.cpp b/dpnp/backend/extensions/lapack/lapack_py.cpp index 4273c90697a4..c0765be7509d 100644 --- a/dpnp/backend/extensions/lapack/lapack_py.cpp +++ b/dpnp/backend/extensions/lapack/lapack_py.cpp @@ -41,8 +41,8 @@ namespace py = pybind11; // populate dispatch vectors void init_dispatch_vectors(void) { - lapack_ext::init_syevd_dispatch_vector(); lapack_ext::init_gesv_dispatch_vector(); + lapack_ext::init_syevd_dispatch_vector(); } // populate dispatch tables diff --git a/dpnp/linalg/dpnp_utils_linalg.py b/dpnp/linalg/dpnp_utils_linalg.py index 7739223b8646..d106fcd6ba9f 100644 --- a/dpnp/linalg/dpnp_utils_linalg.py +++ b/dpnp/linalg/dpnp_utils_linalg.py @@ -51,7 +51,6 @@ def _common_type(*arrays): Key differences from `numpy.common_type`: - It accepts ``bool_`` arrays. - - It doesn't consider ``float16`` arrays since they are not supported. - The default floating-point data type is determined by the capabilities of the device, as indicated by `dpnp.default_float_type()`. @@ -235,7 +234,7 @@ def dpnp_solve(a, b): res_type = _common_type(a, b) if b.size == 0: - return dpnp.empty(b_shape, dtype=res_type) + return dpnp.empty_like(b, dtype=res_type) if a.ndim > 2: reshape = False diff --git a/tests/test_linalg.py b/tests/test_linalg.py index 8162bd4081a8..6230226a0686 100644 --- a/tests/test_linalg.py +++ b/tests/test_linalg.py @@ -5,7 +5,7 @@ import dpnp as inp -from .helper import get_all_dtypes, has_support_aspect64 +from .helper import assert_dtype_allclose, get_all_dtypes, has_support_aspect64 def vvsort(val, vec, size, xp): @@ -534,11 +534,7 @@ def test_solve_diff_type(self, a_dtype, b_dtype): result = inp.linalg.solve(a_dp, b_dp) assert_allclose(expected, result, rtol=1e-06) - - if has_support_aspect64(): - assert expected.dtype == result.dtype - else: - assert expected.dtype.kind == result.dtype.kind + assert_dtype_allclose(result, expected) def test_solve_strides(self): a_np = numpy.array( From cd21c7ffcc5df967747a692cf69ce32a35d36161 Mon Sep 17 00:00:00 2001 From: Vladislav Perevezentsev Date: Tue, 21 Nov 2023 13:40:38 +0100 Subject: [PATCH 33/57] Add res_usm_type variavble and new tests in test_usm_type for dpnp.linalg.solve --- dpnp/linalg/dpnp_utils_linalg.py | 22 +++++++++++++----- tests/test_usm_type.py | 38 ++++++++++++++++++++++++++++++++ 2 files changed, 54 insertions(+), 6 deletions(-) diff --git a/dpnp/linalg/dpnp_utils_linalg.py b/dpnp/linalg/dpnp_utils_linalg.py index d106fcd6ba9f..42e4d33407f0 100644 --- a/dpnp/linalg/dpnp_utils_linalg.py +++ b/dpnp/linalg/dpnp_utils_linalg.py @@ -29,6 +29,7 @@ import dpctl import dpctl.tensor._tensor_impl as ti +import dpctl.utils as du from numpy import prod import dpnp @@ -233,8 +234,9 @@ def dpnp_solve(a, b): ) res_type = _common_type(a, b) + res_usm_type = du.get_coerced_usm_type([a.usm_type, b.usm_type]) if b.size == 0: - return dpnp.empty_like(b, dtype=res_type) + return dpnp.empty_like(b, dtype=res_type, usm_type=res_usm_type) if a.ndim > 2: reshape = False @@ -264,8 +266,12 @@ def dpnp_solve(a, b): # oneMKL LAPACK assumes fortran-like array as input, so # allocate a memory with 'F' order for dpnp array of coefficient matrix # and multiple dependent variables array - coeff_vecs[i] = dpnp.empty_like(a[i], order="F", dtype=res_type) - val_vecs[i] = dpnp.empty_like(b[i], order="F", dtype=res_type) + coeff_vecs[i] = dpnp.empty_like( + a[i], order="F", dtype=res_type, usm_type=res_usm_type + ) + val_vecs[i] = dpnp.empty_like( + b[i], order="F", dtype=res_type, usm_type=res_usm_type + ) # use DPCTL tensor function to fill the coefficient matrix array # and the array of multiple dependent variables with content @@ -297,7 +303,7 @@ def dpnp_solve(a, b): a_ht_copy_ev[i].wait() # combine the list of solutions into a single array - out_v = dpnp.array(val_vecs, order=b_order) + out_v = dpnp.array(val_vecs, order=b_order, usm_type=res_usm_type) if reshape: # shape of the out_v must be equal to the shape of the array of # dependent variables @@ -307,8 +313,12 @@ def dpnp_solve(a, b): # oneMKL LAPACK assumes fortran-like array as input, so # allocate a memory with 'F' order for dpnp array of coefficient matrix # and multiple dependent variables - a_f = dpnp.empty_like(a, order="F", dtype=res_type) - b_f = dpnp.empty_like(b, order="F", dtype=res_type) + a_f = dpnp.empty_like( + a, order="F", dtype=res_type, usm_type=res_usm_type + ) + b_f = dpnp.empty_like( + b, order="F", dtype=res_type, usm_type=res_usm_type + ) # use DPCTL tensor function to fill the coefficient matrix array # and the array of multiple dependent variables with content diff --git a/tests/test_usm_type.py b/tests/test_usm_type.py index 99a39acae887..b60c84e4e50f 100644 --- a/tests/test_usm_type.py +++ b/tests/test_usm_type.py @@ -480,3 +480,41 @@ def test_take(usm_type_x, usm_type_ind): assert x.usm_type == usm_type_x assert ind.usm_type == usm_type_ind assert z.usm_type == du.get_coerced_usm_type([usm_type_x, usm_type_ind]) + + +@pytest.mark.parametrize( + "usm_type_matrix", list_of_usm_types, ids=list_of_usm_types +) +@pytest.mark.parametrize( + "usm_type_vector", list_of_usm_types, ids=list_of_usm_types +) +@pytest.mark.parametrize( + "matrix, vector", + [ + ([[1, 2], [3, 5]], dp.empty((2, 0))), + ([[1, 2], [3, 5]], [1, 2]), + ( + [ + [[1, 1, 1], [0, 2, 5], [2, 5, -1]], + [[3, -1, 1], [1, 2, 3], [2, 3, 1]], + [[1, 4, 1], [1, 2, -2], [4, 1, 2]], + ], + [[6, -4, 27], [9, -6, 15], [15, 1, 11]], + ), + ], + ids=[ + "2D_Matrix_Empty_Vector", + "2D_Matrix_1D_Vector", + "3D_Matrix_and_Vectors", + ], +) +def test_solve(matrix, vector, usm_type_matrix, usm_type_vector): + x = dp.array(matrix, usm_type=usm_type_matrix) + y = dp.array(vector, usm_type=usm_type_vector) + z = dp.linalg.solve(x, y) + + assert x.usm_type == usm_type_matrix + assert y.usm_type == usm_type_vector + assert z.usm_type == du.get_coerced_usm_type( + [usm_type_matrix, usm_type_vector] + ) From 5781f0ce33cfd5366c2f6c06d9fc76eb7636edf1 Mon Sep 17 00:00:00 2001 From: Vladislav Perevezentsev Date: Wed, 22 Nov 2023 13:19:45 +0100 Subject: [PATCH 34/57] Add skipif for test_solve_singular_matrix on cpu --- tests/test_linalg.py | 9 ++++++++- 1 file changed, 8 insertions(+), 1 deletion(-) diff --git a/tests/test_linalg.py b/tests/test_linalg.py index 6230226a0686..ec5dc4287860 100644 --- a/tests/test_linalg.py +++ b/tests/test_linalg.py @@ -5,7 +5,12 @@ import dpnp as inp -from .helper import assert_dtype_allclose, get_all_dtypes, has_support_aspect64 +from .helper import ( + assert_dtype_allclose, + get_all_dtypes, + has_support_aspect64, + is_cpu_device, +) def vvsort(val, vec, size, xp): @@ -561,6 +566,8 @@ def test_solve_strides(self): result = inp.linalg.solve(a_dp[::-2, ::-2], b_dp[::-2]) assert_allclose(expected, result, rtol=1e-05) + # TODO: remove skipif when MKLD-16626 is resolved + @pytest.mark.skipif(is_cpu_device(), reason="MKL bug MKLD-16626") @pytest.mark.parametrize( "matrix, vector", [ From cec81542e448b53843de97b325a6876ede812187 Mon Sep 17 00:00:00 2001 From: Vladislav Perevezentsev Date: Thu, 30 Nov 2023 16:18:05 +0100 Subject: [PATCH 35/57] Modify _common_inexact_type and add a description for it --- dpnp/linalg/dpnp_utils_linalg.py | 31 +++++++++++++++++++++++++------ 1 file changed, 25 insertions(+), 6 deletions(-) diff --git a/dpnp/linalg/dpnp_utils_linalg.py b/dpnp/linalg/dpnp_utils_linalg.py index 42e4d33407f0..7996f3c15242 100644 --- a/dpnp/linalg/dpnp_utils_linalg.py +++ b/dpnp/linalg/dpnp_utils_linalg.py @@ -30,7 +30,7 @@ import dpctl import dpctl.tensor._tensor_impl as ti import dpctl.utils as du -from numpy import prod +from numpy import issubdtype, prod import dpnp import dpnp.backend.extensions.lapack._lapack_impl as li @@ -45,7 +45,7 @@ def _common_type(*arrays): """ _common_type(*arrays) - Common type for linalg + Common type for linear algebra operations. This function determines the common data type for linalg operations. It's designed to be similar in logic to `numpy.linalg.linalg._commonType`. @@ -68,15 +68,34 @@ def _common_type(*arrays): dtypes = [arr.dtype for arr in arrays] - default = dpnp.default_float_type().name - dtype_common = _common_type_internal(default, *dtypes) + default = dpnp.default_float_type() + dtype_common = _common_inexact_type(default, *dtypes) return dtype_common -def _common_type_internal(default_dtype, *dtypes): +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 + 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, + ensuring computational compatibility. + + """ inexact_dtypes = [ - dtype if dtype.kind in "fc" else default_dtype for dtype in dtypes + dt if issubdtype(dt, dpnp.inexact) else default_dtype for dt in dtypes ] return dpnp.result_type(*inexact_dtypes) From c3e5a0fb9002e1fabe6421c8f4bfc00e8d2239a7 Mon Sep 17 00:00:00 2001 From: Vladislav Perevezentsev Date: Thu, 30 Nov 2023 16:31:03 +0100 Subject: [PATCH 36/57] A small update of the desctiption of dpnp.linalg.solve() func --- dpnp/linalg/dpnp_iface_linalg.py | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/dpnp/linalg/dpnp_iface_linalg.py b/dpnp/linalg/dpnp_iface_linalg.py index 5162a4545303..e61194ea84e0 100644 --- a/dpnp/linalg/dpnp_iface_linalg.py +++ b/dpnp/linalg/dpnp_iface_linalg.py @@ -516,7 +516,8 @@ def solve(a, b): Limitations ----------- - Parameter `a` is supported as :class:`dpnp.ndarray` or :class:`dpctl.tensor.usm_ndarray`. + Parameters `a` and `b` are supported as either :class:`dpnp.ndarray` + or :class:`dpctl.tensor.usm_ndarray`. Input array data types are limited by supported DPNP :ref:`Data types`. See Also @@ -532,6 +533,11 @@ def solve(a, b): >>> x array([-1., 1.]) + Check that the solution is correct: + + >>> dp.allclose(dp.dot(a, x), b) + array([ True]) + """ if not dpnp.is_supported_array_type(a): From eb2dd4ced1bacd2755cab48e69ba9617d6ae18ac Mon Sep 17 00:00:00 2001 From: Vladislav Perevezentsev Date: Fri, 1 Dec 2023 16:01:24 +0100 Subject: [PATCH 37/57] Use device param for default_float_type in _common_type --- dpnp/linalg/dpnp_utils_linalg.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/dpnp/linalg/dpnp_utils_linalg.py b/dpnp/linalg/dpnp_utils_linalg.py index 7996f3c15242..3df861149381 100644 --- a/dpnp/linalg/dpnp_utils_linalg.py +++ b/dpnp/linalg/dpnp_utils_linalg.py @@ -52,8 +52,8 @@ def _common_type(*arrays): Key differences from `numpy.common_type`: - It accepts ``bool_`` arrays. - - The default floating-point data type is determined by the capabilities of the device, - as indicated by `dpnp.default_float_type()`. + - 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. @@ -68,7 +68,7 @@ def _common_type(*arrays): dtypes = [arr.dtype for arr in arrays] - default = dpnp.default_float_type() + default = dpnp.default_float_type(device=arrays[0].device) dtype_common = _common_inexact_type(default, *dtypes) return dtype_common From 4780597cb6822ee25d132c0c21ed7aaab44cd80b Mon Sep 17 00:00:00 2001 From: Vladislav Perevezentsev Date: Mon, 4 Dec 2023 11:44:22 +0100 Subject: [PATCH 38/57] Simplify getting 3d array in dpnp_solve --- dpnp/linalg/dpnp_utils_linalg.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/dpnp/linalg/dpnp_utils_linalg.py b/dpnp/linalg/dpnp_utils_linalg.py index 3df861149381..4c57120a9d3f 100644 --- a/dpnp/linalg/dpnp_utils_linalg.py +++ b/dpnp/linalg/dpnp_utils_linalg.py @@ -30,7 +30,7 @@ import dpctl import dpctl.tensor._tensor_impl as ti import dpctl.utils as du -from numpy import issubdtype, prod +from numpy import issubdtype import dpnp import dpnp.backend.extensions.lapack._lapack_impl as li @@ -263,11 +263,11 @@ def dpnp_solve(a, b): if a.ndim > 3: # get 3d input arrays by reshape if a.ndim == b.ndim: - b = b.reshape(prod(b_shape[:-2]), b_shape[-2], b_shape[-1]) + b = b.reshape(-1, b_shape[-2], b_shape[-1]) else: - b = b.reshape(prod(b_shape[:-1]), b_shape[-1]) + b = b.reshape(-1, b_shape[-1]) - a = a.reshape(prod(a_shape[:-2]), a_shape[-2], a_shape[-1]) + a = a.reshape(-1, a_shape[-2], a_shape[-1]) a_usm_arr = dpnp.get_usm_ndarray(a) b_usm_arr = dpnp.get_usm_ndarray(b) From f1b6a81edc5df9841c04ada6e7d98d219f9e202b Mon Sep 17 00:00:00 2001 From: Vladislav Perevezentsev Date: Mon, 4 Dec 2023 12:47:21 +0100 Subject: [PATCH 39/57] Remove unnecessary copying to F order after invoking gesv --- dpnp/linalg/dpnp_utils_linalg.py | 18 ++---------------- 1 file changed, 2 insertions(+), 16 deletions(-) diff --git a/dpnp/linalg/dpnp_utils_linalg.py b/dpnp/linalg/dpnp_utils_linalg.py index 4c57120a9d3f..16c02b1109ed 100644 --- a/dpnp/linalg/dpnp_utils_linalg.py +++ b/dpnp/linalg/dpnp_utils_linalg.py @@ -351,26 +351,12 @@ def dpnp_solve(a, b): # Call the LAPACK extension function _gesv to solve the system of linear # equations with the coefficient square matrix and the dependent variables array. - ht_lapack_ev, lapack_ev = li._gesv( + ht_lapack_ev, _ = li._gesv( exec_q, a_f.get_array(), b_f.get_array(), [a_copy_ev, b_copy_ev] ) - if b_order != "F": - # need to align order of the result of solutions with the - # input array of multiple dependent variables - out_v = dpnp.empty_like(b_f, order=b_order) - ht_copy_out_ev, _ = ti._copy_usm_ndarray_into_usm_ndarray( - src=b_f.get_array(), - dst=out_v.get_array(), - sycl_queue=exec_q, - depends=[lapack_ev], - ) - ht_copy_out_ev.wait() - else: - out_v = b_f - ht_lapack_ev.wait() b_ht_copy_ev.wait() a_ht_copy_ev.wait() - return out_v + return b_f From b8f4cb99f6f809ee6046a0fea4a4b3fa3b500335 Mon Sep 17 00:00:00 2001 From: Vladislav Perevezentsev Date: Mon, 4 Dec 2023 13:03:10 +0100 Subject: [PATCH 40/57] Use get_usm_allocations instead of get_execution_queue --- dpnp/linalg/dpnp_utils_linalg.py | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/dpnp/linalg/dpnp_utils_linalg.py b/dpnp/linalg/dpnp_utils_linalg.py index 16c02b1109ed..76de967288bf 100644 --- a/dpnp/linalg/dpnp_utils_linalg.py +++ b/dpnp/linalg/dpnp_utils_linalg.py @@ -27,13 +27,12 @@ # ***************************************************************************** -import dpctl import dpctl.tensor._tensor_impl as ti -import dpctl.utils as du from numpy import issubdtype import dpnp import dpnp.backend.extensions.lapack._lapack_impl as li +from dpnp.dpnp_utils import get_usm_allocations __all__ = ["dpnp_eigh", "dpnp_solve"] @@ -245,7 +244,7 @@ def dpnp_solve(a, b): a_shape = a.shape b_shape = b.shape - exec_q = dpctl.utils.get_execution_queue((a.sycl_queue, b.sycl_queue)) + res_usm_type, exec_q = get_usm_allocations([a, b]) if exec_q is None: raise ValueError( "Execution placement can not be unambiguously inferred " @@ -253,7 +252,6 @@ def dpnp_solve(a, b): ) res_type = _common_type(a, b) - res_usm_type = du.get_coerced_usm_type([a.usm_type, b.usm_type]) if b.size == 0: return dpnp.empty_like(b, dtype=res_type, usm_type=res_usm_type) From 00ebef19b02c9bb6c324cc5e4cf81af02a8f58ea Mon Sep 17 00:00:00 2001 From: Vladislav Perevezentsev Date: Mon, 4 Dec 2023 13:15:44 +0100 Subject: [PATCH 41/57] Move copying just after the memory allocation --- dpnp/linalg/dpnp_utils_linalg.py | 18 ++++++++++-------- 1 file changed, 10 insertions(+), 8 deletions(-) diff --git a/dpnp/linalg/dpnp_utils_linalg.py b/dpnp/linalg/dpnp_utils_linalg.py index 76de967288bf..5aab461625ad 100644 --- a/dpnp/linalg/dpnp_utils_linalg.py +++ b/dpnp/linalg/dpnp_utils_linalg.py @@ -327,22 +327,24 @@ def dpnp_solve(a, b): out_v = out_v.reshape(orig_shape_b) return out_v else: - # oneMKL LAPACK assumes fortran-like array as input, so - # allocate a memory with 'F' order for dpnp array of coefficient matrix - # and multiple dependent variables + # oneMKL LAPACK gesv overwrites `a` and `b` and assumes fortran-like array as input. + # Allocate 'F' order memory for dpnp arrays to comply with these requirements. a_f = dpnp.empty_like( a, order="F", dtype=res_type, usm_type=res_usm_type ) - b_f = dpnp.empty_like( - b, order="F", dtype=res_type, usm_type=res_usm_type - ) # use DPCTL tensor function to fill the coefficient matrix array - # and the array of multiple dependent variables with content - # from the input arrays + # with content from the input array a a_ht_copy_ev, a_copy_ev = ti._copy_usm_ndarray_into_usm_ndarray( src=a_usm_arr, dst=a_f.get_array(), sycl_queue=a.sycl_queue ) + + b_f = dpnp.empty_like( + b, order="F", dtype=res_type, usm_type=res_usm_type + ) + + # use DPCTL tensor function to fill the array of multiple dependent variables + # with content from the input array b b_ht_copy_ev, b_copy_ev = ti._copy_usm_ndarray_into_usm_ndarray( src=b_usm_arr, dst=b_f.get_array(), sycl_queue=b.sycl_queue ) From 22d4d6f913c53916df6003a0aeb471aa00b90b76 Mon Sep 17 00:00:00 2001 From: Vladislav Perevezentsev Date: Mon, 4 Dec 2023 13:42:01 +0100 Subject: [PATCH 42/57] Add additional checks to gesv implementation --- dpnp/backend/extensions/lapack/gesv.cpp | 13 +++++++++++++ 1 file changed, 13 insertions(+) diff --git a/dpnp/backend/extensions/lapack/gesv.cpp b/dpnp/backend/extensions/lapack/gesv.cpp index f536564fb92f..ffd93130ec0a 100644 --- a/dpnp/backend/extensions/lapack/gesv.cpp +++ b/dpnp/backend/extensions/lapack/gesv.cpp @@ -186,6 +186,13 @@ std::pair ", but a 2-dimensional array is expected."); } + if (dependent_vals_nd != 1 && dependent_vals_nd != 2) { + throw py::value_error( + "The dependent values array has ndim=" + + std::to_string(dependent_vals_nd) + + ", but a 1-dimensional or a 2-dimensional array is expected."); + } + const py::ssize_t *coeff_matrix_shape = coeff_matrix.get_shape_raw(); const py::ssize_t *dependent_vals_shape = dependent_vals.get_shape_raw(); @@ -217,6 +224,12 @@ std::pair "must be F-contiguous"); } + bool is_dependent_vals_f_contig = dependent_vals.is_f_contiguous(); + if (!is_dependent_vals_f_contig) { + throw py::value_error("The array of dependent variables " + "must be F-contiguous"); + } + auto array_types = dpctl_td_ns::usm_ndarray_types(); int coeff_matrix_type_id = array_types.typenum_to_lookup_id(coeff_matrix.get_typenum()); From 04f8f4151d65ae7ef90c9b493692b78ce2772cb9 Mon Sep 17 00:00:00 2001 From: Vladislav Perevezentsev Date: Mon, 4 Dec 2023 16:15:00 +0100 Subject: [PATCH 43/57] Add validation functions for array types and dimensions for linalg funcs --- dpnp/linalg/dpnp_iface_linalg.py | 32 ++++++------------- dpnp/linalg/dpnp_utils_linalg.py | 53 ++++++++++++++++++++++++++++++++ 2 files changed, 63 insertions(+), 22 deletions(-) diff --git a/dpnp/linalg/dpnp_iface_linalg.py b/dpnp/linalg/dpnp_iface_linalg.py index e61194ea84e0..a23dbc8c36e2 100644 --- a/dpnp/linalg/dpnp_iface_linalg.py +++ b/dpnp/linalg/dpnp_iface_linalg.py @@ -47,7 +47,13 @@ from dpnp.dpnp_utils import * from dpnp.linalg.dpnp_algo_linalg import * -from .dpnp_utils_linalg import dpnp_eigh, dpnp_solve +from .dpnp_utils_linalg import ( + _assert_stacked_2d, + _assert_stacked_square, + _assert_supported_array_type, + dpnp_eigh, + dpnp_solve, +) __all__ = [ "cholesky", @@ -540,27 +546,9 @@ def solve(a, b): """ - if not dpnp.is_supported_array_type(a): - raise TypeError( - "An array must be any of supported type, but got {}".format(type(a)) - ) - - if not dpnp.is_supported_array_type(b): - raise TypeError( - "An array must be any of supported type, but got {}".format(type(b)) - ) - - if a.ndim < 2: - raise dpnp.linalg.LinAlgError( - f"{a.ndim}-dimensional array given. The input coefficient " - "array must be at least two-dimensional" - ) - - m, n = a.shape[-2:] - if m != n: - raise dpnp.linalg.LinAlgError( - "Last 2 dimensions of the input coefficient array must be square" - ) + _assert_supported_array_type(a, b) + _assert_stacked_2d(a) + _assert_stacked_square(a) if not ( (a.ndim == b.ndim or a.ndim == b.ndim + 1) diff --git a/dpnp/linalg/dpnp_utils_linalg.py b/dpnp/linalg/dpnp_utils_linalg.py index 5aab461625ad..e865ca134895 100644 --- a/dpnp/linalg/dpnp_utils_linalg.py +++ b/dpnp/linalg/dpnp_utils_linalg.py @@ -40,6 +40,59 @@ _upper_lower = {"U": 0, "L": 1} +def _assert_supported_array_type(*arrays): + """ + Asserts that each array in `arrays` is of a type supported by dpnp. + + Raises `TypeError` if it`s not. + + """ + for a in arrays: + if not dpnp.is_supported_array_type(a): + raise TypeError( + f"The input array must be any of supported type, but got {type(a)}" + ) + + +def _assert_stacked_2d(*arrays): + """ + Asserts that each array in `arrays` is at least two-dimensional. + + Raises `dpnp.linalg.LinAlgError` if it's not. + + """ + for a in arrays: + if a.ndim < 2: + raise dpnp.linalg.LinAlgError( + f"{a.ndim}-dimensional array given. The input " + "array must be at least two-dimensional" + ) + + +def _assert_stacked_square(*arrays): + """ + Asserts that each array in `arrays` is a square matrix. + + Raises `dpnp.linalg.LinAlgError` if any array is not square. + + Precondition: `arrays` are at least 2d. The caller should assert it + beforehand. For example, + + >>> def solve(a): + ... _assert_stacked_2d(a) + ... _assert_stacked_square(a) + ... ... + + """ + + for a in arrays: + m, n = a.shape[-2:] + if m != n: + raise dpnp.linalg.LinAlgError( + "Last 2 dimensions of the input array must be square" + ) + + def _common_type(*arrays): """ _common_type(*arrays) From 366be7a5ad293bceaea58d14da69331ffd08fdea Mon Sep 17 00:00:00 2001 From: Vladislav Perevezentsev Date: Mon, 4 Dec 2023 16:20:16 +0100 Subject: [PATCH 44/57] Update test_solve_diff_type in test_linalg.py --- tests/test_linalg.py | 1 - 1 file changed, 1 deletion(-) diff --git a/tests/test_linalg.py b/tests/test_linalg.py index ec5dc4287860..e0b4f7122b07 100644 --- a/tests/test_linalg.py +++ b/tests/test_linalg.py @@ -538,7 +538,6 @@ def test_solve_diff_type(self, a_dtype, b_dtype): expected = numpy.linalg.solve(a_np, b_np) result = inp.linalg.solve(a_dp, b_dp) - assert_allclose(expected, result, rtol=1e-06) assert_dtype_allclose(result, expected) def test_solve_strides(self): From 08ac7fecb38b6c9fc841542846d9cf08a81f37d0 Mon Sep 17 00:00:00 2001 From: Vladislav Perevezentsev Date: Tue, 5 Dec 2023 15:18:23 +0100 Subject: [PATCH 45/57] Address the remarks --- dpnp/backend/extensions/lapack/gesv.cpp | 2 +- dpnp/linalg/dpnp_utils_linalg.py | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/dpnp/backend/extensions/lapack/gesv.cpp b/dpnp/backend/extensions/lapack/gesv.cpp index ffd93130ec0a..72e5aa806714 100644 --- a/dpnp/backend/extensions/lapack/gesv.cpp +++ b/dpnp/backend/extensions/lapack/gesv.cpp @@ -186,7 +186,7 @@ std::pair ", but a 2-dimensional array is expected."); } - if (dependent_vals_nd != 1 && dependent_vals_nd != 2) { + if (dependent_vals_nd > 2) { throw py::value_error( "The dependent values array has ndim=" + std::to_string(dependent_vals_nd) + diff --git a/dpnp/linalg/dpnp_utils_linalg.py b/dpnp/linalg/dpnp_utils_linalg.py index e865ca134895..6a91bd1cd584 100644 --- a/dpnp/linalg/dpnp_utils_linalg.py +++ b/dpnp/linalg/dpnp_utils_linalg.py @@ -387,7 +387,7 @@ def dpnp_solve(a, b): ) # use DPCTL tensor function to fill the coefficient matrix array - # with content from the input array a + # with content from the input array `a`` a_ht_copy_ev, a_copy_ev = ti._copy_usm_ndarray_into_usm_ndarray( src=a_usm_arr, dst=a_f.get_array(), sycl_queue=a.sycl_queue ) @@ -397,7 +397,7 @@ def dpnp_solve(a, b): ) # use DPCTL tensor function to fill the array of multiple dependent variables - # with content from the input array b + # with content from the input array `b` b_ht_copy_ev, b_copy_ev = ti._copy_usm_ndarray_into_usm_ndarray( src=b_usm_arr, dst=b_f.get_array(), sycl_queue=b.sycl_queue ) From eb6a8408c57f7668538119e00a82556d5ab77244 Mon Sep 17 00:00:00 2001 From: Vladislav Perevezentsev Date: Tue, 5 Dec 2023 17:55:04 +0100 Subject: [PATCH 46/57] Small update --- dpnp/linalg/dpnp_utils_linalg.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/dpnp/linalg/dpnp_utils_linalg.py b/dpnp/linalg/dpnp_utils_linalg.py index 6a91bd1cd584..571422e04120 100644 --- a/dpnp/linalg/dpnp_utils_linalg.py +++ b/dpnp/linalg/dpnp_utils_linalg.py @@ -387,7 +387,7 @@ def dpnp_solve(a, b): ) # use DPCTL tensor function to fill the coefficient matrix array - # with content from the input array `a`` + # with content from the input array `a` a_ht_copy_ev, a_copy_ev = ti._copy_usm_ndarray_into_usm_ndarray( src=a_usm_arr, dst=a_f.get_array(), sycl_queue=a.sycl_queue ) From 61a60733aa2ba029e5bd98dbe25b79e0b394f1ab Mon Sep 17 00:00:00 2001 From: Vladislav Perevezentsev Date: Tue, 5 Dec 2023 17:55:26 +0100 Subject: [PATCH 47/57] qwe --- qwe.py | 16 ++++++++++++++++ tests/helper.py | 17 +++++++++++++++++ tests/test_sum.py | 9 ++++++++- 3 files changed, 41 insertions(+), 1 deletion(-) create mode 100644 qwe.py diff --git a/qwe.py b/qwe.py new file mode 100644 index 000000000000..41bd640e6583 --- /dev/null +++ b/qwe.py @@ -0,0 +1,16 @@ +import dpnp +import numpy +from tests.helper import assert_dtype_allclose + + +dtype='float32' +data = numpy.arange(100, dtype=dtype) +dpnp_data = dpnp.array(data, device='cpu') + +np_res = numpy.fft.fft(data) +dpnp_res = dpnp.fft.fft(dpnp_data) + +print("numpy: ",np_res.dtype) +print("dpnp: ",dpnp_res.dtype) + +assert_dtype_allclose(dpnp_res, np_res) diff --git a/tests/helper.py b/tests/helper.py index b3d816e769ac..04abf80f781a 100644 --- a/tests/helper.py +++ b/tests/helper.py @@ -15,6 +15,7 @@ def assert_dtype_allclose(dpnp_arr, numpy_arr, check_type=True): """ + list_64bit_types = [numpy.float64, numpy.complex128] is_inexact = lambda x: dpnp.issubdtype(x.dtype, dpnp.inexact) if is_inexact(dpnp_arr) or is_inexact(numpy_arr): tol = 8 * max( @@ -22,6 +23,22 @@ def assert_dtype_allclose(dpnp_arr, numpy_arr, check_type=True): numpy.finfo(numpy_arr.dtype).resolution, ) assert_allclose(dpnp_arr.asnumpy(), numpy_arr, atol=tol, rtol=tol) + if check_type: + numpy_arr_dtype = numpy_arr.dtype + dpnp_arr_dtype = dpnp_arr.dtype + dpnp_arr_dev = dpnp_arr.sycl_device + is_np_arr_f2 = numpy_arr_dtype == numpy.float16 + + if is_np_arr_f2: + if has_support_aspect16(dpnp_arr_dev): + assert dpnp_arr_dtype == numpy_arr_dtype + elif ( + numpy_arr_dtype not in list_64bit_types + or has_support_aspect64(dpnp_arr_dev) + ): + assert dpnp_arr_dtype == numpy_arr_dtype + else: + assert dpnp_arr_dtype.kind == numpy_arr_dtype.kind else: assert_array_equal(dpnp_arr.asnumpy(), numpy_arr) if check_type: diff --git a/tests/test_sum.py b/tests/test_sum.py index 16b17847f270..f34d1cf6252d 100644 --- a/tests/test_sum.py +++ b/tests/test_sum.py @@ -27,10 +27,17 @@ def test_sum_float(dtype): ) ia = dpnp.array(a) + # Flag for type check in special cases + # Skip dtype checks when dpnp handles float32 arrays + # as `dpnp.sum()` and `numpy.sum()` return different dtypes + check_dtype = dtype != dpnp.float32 for axis in range(len(a)): result = dpnp.sum(ia, axis=axis) expected = numpy.sum(a, axis=axis) - assert_dtype_allclose(result, expected) + assert_dtype_allclose(result, expected, check_type=check_dtype) + if not check_dtype: + # Ensure dtype kind matches when check_dtype is False + assert result.dtype.kind == expected.dtype.kind def test_sum_int(): From eba811a0c26f348dc069ecc8737903b9b92486db Mon Sep 17 00:00:00 2001 From: Vladislav Perevezentsev Date: Thu, 7 Dec 2023 15:16:27 +0100 Subject: [PATCH 48/57] Rename assert funcs and make them external in dpnp_utils_linalg --- dpnp/linalg/dpnp_iface_linalg.py | 11 +++--- dpnp/linalg/dpnp_utils_linalg.py | 62 +++++++++++++++++++++----------- 2 files changed, 47 insertions(+), 26 deletions(-) diff --git a/dpnp/linalg/dpnp_iface_linalg.py b/dpnp/linalg/dpnp_iface_linalg.py index 34c33db2bf35..c9c9c855728a 100644 --- a/dpnp/linalg/dpnp_iface_linalg.py +++ b/dpnp/linalg/dpnp_iface_linalg.py @@ -48,9 +48,8 @@ from dpnp.linalg.dpnp_algo_linalg import * from .dpnp_utils_linalg import ( - _assert_stacked_2d, - _assert_stacked_square, - _assert_supported_array_type, + check_stacked_2d, + check_stacked_square, dpnp_eigh, dpnp_solve, ) @@ -543,9 +542,9 @@ def solve(a, b): """ - _assert_supported_array_type(a, b) - _assert_stacked_2d(a) - _assert_stacked_square(a) + dpnp.check_supported_arrays_type(a, b) + check_stacked_2d(a) + check_stacked_square(a) if not ( (a.ndim == b.ndim or a.ndim == b.ndim + 1) diff --git a/dpnp/linalg/dpnp_utils_linalg.py b/dpnp/linalg/dpnp_utils_linalg.py index 571422e04120..aab04f630c7c 100644 --- a/dpnp/linalg/dpnp_utils_linalg.py +++ b/dpnp/linalg/dpnp_utils_linalg.py @@ -34,33 +34,40 @@ import dpnp.backend.extensions.lapack._lapack_impl as li from dpnp.dpnp_utils import get_usm_allocations -__all__ = ["dpnp_eigh", "dpnp_solve"] +__all__ = [ + "check_stacked_2d", + "check_stacked_square", + "dpnp_eigh", + "dpnp_solve", +] _jobz = {"N": 0, "V": 1} _upper_lower = {"U": 0, "L": 1} -def _assert_supported_array_type(*arrays): +def check_stacked_2d(*arrays): """ - Asserts that each array in `arrays` is of a type supported by dpnp. + Return ``True`` if each array in `arrays` has at least two dimensions. - Raises `TypeError` if it`s not. + If any array is less than two-dimensional, `dpnp.linalg.LinAlgError` will be raised. - """ - for a in arrays: - if not dpnp.is_supported_array_type(a): - raise TypeError( - f"The input array must be any of supported type, but got {type(a)}" - ) + Parameters + ---------- + arrays : {dpnp_array, usm_ndarray} + A sequence of input arrays to check for dimensionality. + Returns + ------- + out : bool + ``True`` if each array in `arrays` is at least two-dimensional. -def _assert_stacked_2d(*arrays): - """ - Asserts that each array in `arrays` is at least two-dimensional. - - Raises `dpnp.linalg.LinAlgError` if it's not. + Raises + ------ + dpnp.linalg.LinAlgError + If any array in `arrays` is less than two-dimensional. """ + for a in arrays: if a.ndim < 2: raise dpnp.linalg.LinAlgError( @@ -69,20 +76,35 @@ def _assert_stacked_2d(*arrays): ) -def _assert_stacked_square(*arrays): +def check_stacked_square(*arrays): """ - Asserts that each array in `arrays` is a square matrix. + Return ``True`` if each array in `arrays` is a square matrix. - Raises `dpnp.linalg.LinAlgError` if any array is not square. + If any array does not form a square matrix, `dpnp.linalg.LinAlgError` will be raised. Precondition: `arrays` are at least 2d. The caller should assert it beforehand. For example, >>> def solve(a): - ... _assert_stacked_2d(a) - ... _assert_stacked_square(a) + ... check_stacked_2d(a) + ... check_stacked_square(a) ... ... + Parameters + ---------- + arrays : {dpnp_array, usm_ndarray} + A sequence of input arrays to check for square matrix shape. + + Returns + ------- + out : bool + ``True`` if each array in `arrays` forms a square matrix. + + Raises + ------ + dpnp.linalg.LinAlgError + If any array in `arrays` does not form a square matrix. + """ for a in arrays: From 3a9d459959df7dc3433d60f79fc3c9311a29655c Mon Sep 17 00:00:00 2001 From: Vladislav Perevezentsev Date: Thu, 7 Dec 2023 15:17:08 +0100 Subject: [PATCH 49/57] Use assert_dtype_allclose for test_solve in test_sycl_queue --- tests/test_sycl_queue.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/test_sycl_queue.py b/tests/test_sycl_queue.py index 4a6b839c45fd..95b9d458c523 100644 --- a/tests/test_sycl_queue.py +++ b/tests/test_sycl_queue.py @@ -1338,7 +1338,7 @@ def test_solve(device): result = dpnp.linalg.solve(dpnp_x, dpnp_y) expected = numpy.linalg.solve(numpy_x, numpy_y) - assert_allclose(expected, result, rtol=1e-06) + assert_dtype_allclose(result, expected) result_queue = result.sycl_queue From 0de696897c2becff12a7002ce6aa55b4e0f30ae5 Mon Sep 17 00:00:00 2001 From: Vladislav Perevezentsev Date: Thu, 7 Dec 2023 15:21:28 +0100 Subject: [PATCH 50/57] Remove an unnecessary file --- qwe.py | 16 ---------------- 1 file changed, 16 deletions(-) delete mode 100644 qwe.py diff --git a/qwe.py b/qwe.py deleted file mode 100644 index 41bd640e6583..000000000000 --- a/qwe.py +++ /dev/null @@ -1,16 +0,0 @@ -import dpnp -import numpy -from tests.helper import assert_dtype_allclose - - -dtype='float32' -data = numpy.arange(100, dtype=dtype) -dpnp_data = dpnp.array(data, device='cpu') - -np_res = numpy.fft.fft(data) -dpnp_res = dpnp.fft.fft(dpnp_data) - -print("numpy: ",np_res.dtype) -print("dpnp: ",dpnp_res.dtype) - -assert_dtype_allclose(dpnp_res, np_res) From df72c77b781ad859b15681a0031ebd0f0682bc48 Mon Sep 17 00:00:00 2001 From: Vladislav Perevezentsev Date: Thu, 7 Dec 2023 18:09:08 +0100 Subject: [PATCH 51/57] Fix validation for CI --- tests/test_linalg.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/test_linalg.py b/tests/test_linalg.py index e0b4f7122b07..6c8c2a5f04e3 100644 --- a/tests/test_linalg.py +++ b/tests/test_linalg.py @@ -566,7 +566,7 @@ def test_solve_strides(self): assert_allclose(expected, result, rtol=1e-05) # TODO: remove skipif when MKLD-16626 is resolved - @pytest.mark.skipif(is_cpu_device(), reason="MKL bug MKLD-16626") + @pytest.mark.skipif(is_cpu_device(), reason="MKLD-16626") @pytest.mark.parametrize( "matrix, vector", [ From 4991a819088455a52909265b6ff8b5cbe7fa38f2 Mon Sep 17 00:00:00 2001 From: Vladislav Perevezentsev Date: Fri, 8 Dec 2023 11:48:05 +0100 Subject: [PATCH 52/57] Remove eqec_q check that will never happen --- dpnp/linalg/dpnp_utils_linalg.py | 5 ----- 1 file changed, 5 deletions(-) diff --git a/dpnp/linalg/dpnp_utils_linalg.py b/dpnp/linalg/dpnp_utils_linalg.py index aab04f630c7c..fd6353358c28 100644 --- a/dpnp/linalg/dpnp_utils_linalg.py +++ b/dpnp/linalg/dpnp_utils_linalg.py @@ -320,11 +320,6 @@ def dpnp_solve(a, b): b_shape = b.shape res_usm_type, exec_q = get_usm_allocations([a, b]) - if exec_q is None: - raise ValueError( - "Execution placement can not be unambiguously inferred " - "from input arguments." - ) res_type = _common_type(a, b) if b.size == 0: From 9b5a5a58dea8a90b1d106c0b7e8978c96413abed Mon Sep 17 00:00:00 2001 From: Vladislav Perevezentsev Date: Mon, 11 Dec 2023 11:34:21 +0100 Subject: [PATCH 53/57] Set usm_type for out_v --- dpnp/linalg/dpnp_utils_linalg.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/dpnp/linalg/dpnp_utils_linalg.py b/dpnp/linalg/dpnp_utils_linalg.py index fd6353358c28..a3a2802072c9 100644 --- a/dpnp/linalg/dpnp_utils_linalg.py +++ b/dpnp/linalg/dpnp_utils_linalg.py @@ -390,7 +390,9 @@ def dpnp_solve(a, b): a_ht_copy_ev[i].wait() # combine the list of solutions into a single array - out_v = dpnp.array(val_vecs, order=b_order, usm_type=res_usm_type) + out_v = dpnp.array( + val_vecs, order=b_order, dtype=res_type, usm_type=res_usm_type + ) if reshape: # shape of the out_v must be equal to the shape of the array of # dependent variables From 3c7ad077795b7b9a7b1d1d857860c2925a2f9559 Mon Sep 17 00:00:00 2001 From: Vladislav Perevezentsev Date: Mon, 11 Dec 2023 11:53:30 +0100 Subject: [PATCH 54/57] Update test_solve_singular_empty --- tests/third_party/cupy/linalg_tests/test_solve.py | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/tests/third_party/cupy/linalg_tests/test_solve.py b/tests/third_party/cupy/linalg_tests/test_solve.py index 5de4aacc789b..718fd03400ec 100644 --- a/tests/third_party/cupy/linalg_tests/test_solve.py +++ b/tests/third_party/cupy/linalg_tests/test_solve.py @@ -65,9 +65,10 @@ def test_solve_singular_empty(self): for xp in (numpy, cupy): a = xp.zeros((3, 3)) # singular b = xp.empty((3, 0)) # nrhs = 0 - # numpy <= 1.24.* raises LinAlgError when b.size == 0 - # numpy >= 1.25 returns an empty array - if xp == numpy and testing.numpy_satisfies("<1.25.0"): + # Different behavior of numpy and dpnp in this case + # numpy raises LinAlgError when b.size == 0 and a is a singular matrix + # dpnp returns an empty array + if xp == numpy: with pytest.raises(numpy.linalg.LinAlgError): xp.linalg.solve(a, b) else: From 3001872e3e12047c304b19fbe35db0e08579153f Mon Sep 17 00:00:00 2001 From: Vladislav Perevezentsev Date: Mon, 11 Dec 2023 17:58:03 +0100 Subject: [PATCH 55/57] Skip test_solve_singular_empty --- .../cupy/linalg_tests/test_solve.py | 21 +++++++------------ 1 file changed, 8 insertions(+), 13 deletions(-) diff --git a/tests/third_party/cupy/linalg_tests/test_solve.py b/tests/third_party/cupy/linalg_tests/test_solve.py index 718fd03400ec..5054721d5dc7 100644 --- a/tests/third_party/cupy/linalg_tests/test_solve.py +++ b/tests/third_party/cupy/linalg_tests/test_solve.py @@ -61,19 +61,14 @@ def check_shape(self, a_shape, b_shape, error_types): with pytest.raises(error_type): xp.linalg.solve(a, b) - def test_solve_singular_empty(self): - for xp in (numpy, cupy): - a = xp.zeros((3, 3)) # singular - b = xp.empty((3, 0)) # nrhs = 0 - # Different behavior of numpy and dpnp in this case - # numpy raises LinAlgError when b.size == 0 and a is a singular matrix - # dpnp returns an empty array - if xp == numpy: - with pytest.raises(numpy.linalg.LinAlgError): - xp.linalg.solve(a, b) - else: - result = xp.linalg.solve(a, b) - assert result.size == 0 + # Undefined behavior is implementation-dependent: + # Numpy with OpenBLAS returns an empty array while numpy with MKL raises LinAlgError + @pytest.mark.skip("Undefined behavior") + def test_solve_singular_empty(self, xp): + a = xp.zeros((3, 3)) # singular + b = xp.empty((3, 0)) # nrhs = 0 + # LinAlgError("Singular matrix") is not raised + return xp.linalg.solve(a, b) def test_invalid_shape(self): linalg_errors = { From 5392a45092b240724cd481df809f1762c8e00d6f Mon Sep 17 00:00:00 2001 From: Vladislav Perevezentsev Date: Tue, 12 Dec 2023 16:12:53 +0100 Subject: [PATCH 56/57] Fix validation fell --- tests/third_party/cupy/linalg_tests/test_solve.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/tests/third_party/cupy/linalg_tests/test_solve.py b/tests/third_party/cupy/linalg_tests/test_solve.py index 5054721d5dc7..b40f2ed4ea0f 100644 --- a/tests/third_party/cupy/linalg_tests/test_solve.py +++ b/tests/third_party/cupy/linalg_tests/test_solve.py @@ -62,7 +62,8 @@ def check_shape(self, a_shape, b_shape, error_types): xp.linalg.solve(a, b) # Undefined behavior is implementation-dependent: - # Numpy with OpenBLAS returns an empty array while numpy with MKL raises LinAlgError + # Numpy with OpenBLAS returns an empty array + # while numpy with Intel MKL raises LinAlgError @pytest.mark.skip("Undefined behavior") def test_solve_singular_empty(self, xp): a = xp.zeros((3, 3)) # singular From c105570a3ebd335f7969383bfba684f7af686a6a Mon Sep 17 00:00:00 2001 From: Vladislav Perevezentsev Date: Tue, 12 Dec 2023 16:58:23 +0100 Subject: [PATCH 57/57] A small update --- tests/third_party/cupy/linalg_tests/test_solve.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/third_party/cupy/linalg_tests/test_solve.py b/tests/third_party/cupy/linalg_tests/test_solve.py index b40f2ed4ea0f..6194cf6b8ac4 100644 --- a/tests/third_party/cupy/linalg_tests/test_solve.py +++ b/tests/third_party/cupy/linalg_tests/test_solve.py @@ -63,7 +63,7 @@ def check_shape(self, a_shape, b_shape, error_types): # Undefined behavior is implementation-dependent: # Numpy with OpenBLAS returns an empty array - # while numpy with Intel MKL raises LinAlgError + # while numpy with OneMKL raises LinAlgError @pytest.mark.skip("Undefined behavior") def test_solve_singular_empty(self, xp): a = xp.zeros((3, 3)) # singular