diff --git a/dpnp/CMakeLists.txt b/dpnp/CMakeLists.txt index 6ad063a6776..9c79d5af385 100644 --- a/dpnp/CMakeLists.txt +++ b/dpnp/CMakeLists.txt @@ -64,5 +64,4 @@ add_subdirectory(backend/extensions/sycl_ext) add_subdirectory(dpnp_algo) add_subdirectory(dpnp_utils) add_subdirectory(fft) -add_subdirectory(linalg) add_subdirectory(random) diff --git a/dpnp/backend/include/dpnp_iface_fptr.hpp b/dpnp/backend/include/dpnp_iface_fptr.hpp index 4825ef4cd17..b63f0e00431 100644 --- a/dpnp/backend/include/dpnp_iface_fptr.hpp +++ b/dpnp/backend/include/dpnp_iface_fptr.hpp @@ -126,11 +126,7 @@ enum class DPNPFuncName : size_t DPNP_FN_EDIFF1D_EXT, /**< Used in numpy.ediff1d() impl, requires extra parameters */ DPNP_FN_EIG, /**< Used in numpy.linalg.eig() impl */ - DPNP_FN_EIG_EXT, /**< Used in numpy.linalg.eig() impl, requires extra - parameters */ DPNP_FN_EIGVALS, /**< Used in numpy.linalg.eigvals() impl */ - DPNP_FN_EIGVALS_EXT, /**< Used in numpy.linalg.eigvals() impl, requires - extra parameters */ DPNP_FN_ERF, /**< Used in scipy.special.erf impl */ DPNP_FN_ERF_EXT, /**< Used in scipy.special.erf impl, requires extra parameters */ diff --git a/dpnp/backend/kernels/dpnp_krnl_common.cpp b/dpnp/backend/kernels/dpnp_krnl_common.cpp index 04eac54310d..423851e4bfd 100644 --- a/dpnp/backend/kernels/dpnp_krnl_common.cpp +++ b/dpnp/backend/kernels/dpnp_krnl_common.cpp @@ -630,15 +630,6 @@ template void (*dpnp_eig_default_c)(const void *, void *, void *, size_t) = dpnp_eig_c<_DataType, _ResultType>; -template -DPCTLSyclEventRef (*dpnp_eig_ext_c)(DPCTLSyclQueueRef, - const void *, - void *, - void *, - size_t, - const DPCTLEventVectorRef) = - dpnp_eig_c<_DataType, _ResultType>; - template DPCTLSyclEventRef dpnp_eigvals_c(DPCTLSyclQueueRef q_ref, const void *array_in, @@ -724,14 +715,6 @@ void (*dpnp_eigvals_default_c)(const void *, void *, size_t) = dpnp_eigvals_c<_DataType, _ResultType>; -template -DPCTLSyclEventRef (*dpnp_eigvals_ext_c)(DPCTLSyclQueueRef, - const void *, - void *, - size_t, - const DPCTLEventVectorRef) = - dpnp_eigvals_c<_DataType, _ResultType>; - template class dpnp_initval_c_kernel; @@ -1083,27 +1066,6 @@ void func_map_init_linalg(func_map_t &fmap) fmap[DPNPFuncName::DPNP_FN_EIG][eft_DBL][eft_DBL] = { eft_DBL, (void *)dpnp_eig_default_c}; - fmap[DPNPFuncName::DPNP_FN_EIG_EXT][eft_INT][eft_INT] = { - get_default_floating_type(), - (void *)dpnp_eig_ext_c< - int32_t, func_type_map_t::find_type>, - get_default_floating_type(), - (void *)dpnp_eig_ext_c< - int32_t, func_type_map_t::find_type< - get_default_floating_type()>>}; - fmap[DPNPFuncName::DPNP_FN_EIG_EXT][eft_LNG][eft_LNG] = { - get_default_floating_type(), - (void *)dpnp_eig_ext_c< - int64_t, func_type_map_t::find_type>, - get_default_floating_type(), - (void *)dpnp_eig_ext_c< - int64_t, func_type_map_t::find_type< - get_default_floating_type()>>}; - fmap[DPNPFuncName::DPNP_FN_EIG_EXT][eft_FLT][eft_FLT] = { - eft_FLT, (void *)dpnp_eig_ext_c}; - fmap[DPNPFuncName::DPNP_FN_EIG_EXT][eft_DBL][eft_DBL] = { - eft_DBL, (void *)dpnp_eig_ext_c}; - fmap[DPNPFuncName::DPNP_FN_EIGVALS][eft_INT][eft_INT] = { eft_DBL, (void *)dpnp_eigvals_default_c}; fmap[DPNPFuncName::DPNP_FN_EIGVALS][eft_LNG][eft_LNG] = { @@ -1113,27 +1075,6 @@ void func_map_init_linalg(func_map_t &fmap) fmap[DPNPFuncName::DPNP_FN_EIGVALS][eft_DBL][eft_DBL] = { eft_DBL, (void *)dpnp_eigvals_default_c}; - fmap[DPNPFuncName::DPNP_FN_EIGVALS_EXT][eft_INT][eft_INT] = { - get_default_floating_type(), - (void *)dpnp_eigvals_ext_c< - int32_t, func_type_map_t::find_type>, - get_default_floating_type(), - (void *)dpnp_eigvals_ext_c< - int32_t, func_type_map_t::find_type< - get_default_floating_type()>>}; - fmap[DPNPFuncName::DPNP_FN_EIGVALS_EXT][eft_LNG][eft_LNG] = { - get_default_floating_type(), - (void *)dpnp_eigvals_ext_c< - int64_t, func_type_map_t::find_type>, - get_default_floating_type(), - (void *)dpnp_eigvals_ext_c< - int64_t, func_type_map_t::find_type< - get_default_floating_type()>>}; - fmap[DPNPFuncName::DPNP_FN_EIGVALS_EXT][eft_FLT][eft_FLT] = { - eft_FLT, (void *)dpnp_eigvals_ext_c}; - fmap[DPNPFuncName::DPNP_FN_EIGVALS_EXT][eft_DBL][eft_DBL] = { - eft_DBL, (void *)dpnp_eigvals_ext_c}; - fmap[DPNPFuncName::DPNP_FN_INITVAL][eft_BLN][eft_BLN] = { eft_BLN, (void *)dpnp_initval_default_c}; fmap[DPNPFuncName::DPNP_FN_INITVAL][eft_INT][eft_INT] = { diff --git a/dpnp/dpnp_algo/dpnp_algo.pxd b/dpnp/dpnp_algo/dpnp_algo.pxd index ee12ffa4947..5859347aac4 100644 --- a/dpnp/dpnp_algo/dpnp_algo.pxd +++ b/dpnp/dpnp_algo/dpnp_algo.pxd @@ -42,8 +42,6 @@ cdef extern from "dpnp_iface_fptr.hpp" namespace "DPNPFuncName": # need this na DPNP_FN_DIAG_INDICES_EXT DPNP_FN_DIAGONAL_EXT DPNP_FN_EDIFF1D_EXT - DPNP_FN_EIG_EXT - DPNP_FN_EIGVALS_EXT DPNP_FN_ERF_EXT DPNP_FN_FABS_EXT DPNP_FN_FFT_FFT_EXT diff --git a/dpnp/linalg/CMakeLists.txt b/dpnp/linalg/CMakeLists.txt deleted file mode 100644 index a04d5f3b64e..00000000000 --- a/dpnp/linalg/CMakeLists.txt +++ /dev/null @@ -1,7 +0,0 @@ -# Building dpnp_algo_linalg Cython extension - -build_dpnp_cython_ext_with_backend( - dpnp_algo_linalg - ${CMAKE_CURRENT_SOURCE_DIR}/dpnp_algo_linalg.pyx - dpnp/linalg - ) diff --git a/dpnp/linalg/dpnp_algo_linalg.pyx b/dpnp/linalg/dpnp_algo_linalg.pyx deleted file mode 100644 index 367247c2368..00000000000 --- a/dpnp/linalg/dpnp_algo_linalg.pyx +++ /dev/null @@ -1,147 +0,0 @@ -# cython: language_level=3 -# cython: linetrace=True -# -*- coding: utf-8 -*- -# ***************************************************************************** -# Copyright (c) 2016-2024, Intel Corporation -# All rights reserved. -# -# Redistribution and use in source and binary forms, with or without -# modification, are permitted provided that the following conditions are met: -# - Redistributions of source code must retain the above copyright notice, -# this list of conditions and the following disclaimer. -# - Redistributions in binary form must reproduce the above copyright notice, -# this list of conditions and the following disclaimer in the documentation -# and/or other materials provided with the distribution. -# -# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" -# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE -# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE -# ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE -# LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR -# CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF -# SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS -# INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN -# CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) -# ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF -# THE POSSIBILITY OF SUCH DAMAGE. -# ***************************************************************************** - -"""Module Backend - -This module contains interface functions between C backend layer -and the rest of the library - -""" - -import numpy - -from dpnp.dpnp_algo cimport * - -import dpnp -import dpnp.dpnp_utils as utils_py - -cimport numpy - -cimport dpnp.dpnp_utils as utils - -__all__ = [ - "dpnp_eig", - "dpnp_eigvals", -] - - -# C function pointer to the C library template functions -ctypedef c_dpctl.DPCTLSyclEventRef(*custom_linalg_1in_1out_with_size_func_ptr_t_)(c_dpctl.DPCTLSyclQueueRef, - void *, void * , size_t, - const c_dpctl.DPCTLEventVectorRef) -ctypedef c_dpctl.DPCTLSyclEventRef(*custom_linalg_2in_1out_func_ptr_t)(c_dpctl.DPCTLSyclQueueRef, - void *, void * , void * , size_t, - const c_dpctl.DPCTLEventVectorRef) - - -cpdef tuple dpnp_eig(utils.dpnp_descriptor x1): - cdef shape_type_c x1_shape = x1.shape - - cdef size_t size = 0 if x1_shape.empty() else x1_shape.front() - - cdef DPNPFuncType param1_type = dpnp_dtype_to_DPNPFuncType(x1.dtype) - cdef DPNPFuncData kernel_data = get_dpnp_function_ptr(DPNP_FN_EIG_EXT, param1_type, param1_type) - - x1_obj = x1.get_array() - - cdef (DPNPFuncType, void *) ret_type_and_func = utils.get_ret_type_and_func(kernel_data, - x1_obj.sycl_device.has_aspect_fp64) - cdef DPNPFuncType return_type = ret_type_and_func[0] - cdef custom_linalg_2in_1out_func_ptr_t func = < custom_linalg_2in_1out_func_ptr_t > ret_type_and_func[1] - - cdef utils.dpnp_descriptor res_val = utils.create_output_descriptor((size,), - return_type, - None, - device=x1_obj.sycl_device, - usm_type=x1_obj.usm_type, - sycl_queue=x1_obj.sycl_queue) - cdef utils.dpnp_descriptor res_vec = utils.create_output_descriptor(x1_shape, - return_type, - None, - device=x1_obj.sycl_device, - usm_type=x1_obj.usm_type, - sycl_queue=x1_obj.sycl_queue) - - result_sycl_queue = res_val.get_array().sycl_queue - - cdef c_dpctl.SyclQueue q = result_sycl_queue - cdef c_dpctl.DPCTLSyclQueueRef q_ref = q.get_queue_ref() - - # call FPTR function - cdef c_dpctl.DPCTLSyclEventRef event_ref = func(q_ref, - x1.get_data(), - res_val.get_data(), - res_vec.get_data(), - size, - NULL) # dep_events_ref - - with nogil: c_dpctl.DPCTLEvent_WaitAndThrow(event_ref) - c_dpctl.DPCTLEvent_Delete(event_ref) - - return (res_val.get_pyobj(), res_vec.get_pyobj()) - - -cpdef utils.dpnp_descriptor dpnp_eigvals(utils.dpnp_descriptor input): - cdef shape_type_c input_shape = input.shape - - cdef size_t size = 0 if input_shape.empty() else input_shape.front() - - cdef DPNPFuncType param1_type = dpnp_dtype_to_DPNPFuncType(input.dtype) - cdef DPNPFuncData kernel_data = get_dpnp_function_ptr(DPNP_FN_EIGVALS_EXT, param1_type, param1_type) - - input_obj = input.get_array() - - cdef (DPNPFuncType, void *) ret_type_and_func = utils.get_ret_type_and_func(kernel_data, - input_obj.sycl_device.has_aspect_fp64) - cdef DPNPFuncType return_type = ret_type_and_func[0] - cdef custom_linalg_1in_1out_with_size_func_ptr_t_ func = < custom_linalg_1in_1out_with_size_func_ptr_t_ > ret_type_and_func[1] - - # create result array with type given by FPTR data - cdef utils.dpnp_descriptor res_val = utils.create_output_descriptor((size,), - return_type, - None, - device=input_obj.sycl_device, - usm_type=input_obj.usm_type, - sycl_queue=input_obj.sycl_queue) - - result_sycl_queue = res_val.get_array().sycl_queue - - cdef c_dpctl.SyclQueue q = result_sycl_queue - cdef c_dpctl.DPCTLSyclQueueRef q_ref = q.get_queue_ref() - - # call FPTR function - cdef c_dpctl.DPCTLSyclEventRef event_ref = func(q_ref, - input.get_data(), - res_val.get_data(), - size, - NULL) # dep_events_ref - - with nogil: c_dpctl.DPCTLEvent_WaitAndThrow(event_ref) - c_dpctl.DPCTLEvent_Delete(event_ref) - - return res_val diff --git a/dpnp/linalg/dpnp_iface_linalg.py b/dpnp/linalg/dpnp_iface_linalg.py index 27dd7408191..1e7f6cb0bb3 100644 --- a/dpnp/linalg/dpnp_iface_linalg.py +++ b/dpnp/linalg/dpnp_iface_linalg.py @@ -40,9 +40,6 @@ import numpy import dpnp -from dpnp.dpnp_algo import * -from dpnp.dpnp_utils import * -from dpnp.linalg.dpnp_algo_linalg import * from .dpnp_utils_linalg import ( check_stacked_2d, @@ -249,20 +246,104 @@ def det(a): return dpnp_det(a) -def eig(x1): +def eig(a): """ Compute the eigenvalues and right eigenvectors of a square array. - .. seealso:: :obj:`numpy.linalg.eig` + For full documentation refer to :obj:`numpy.linalg.eig`. + + Parameters + ---------- + a : (..., M, M) {dpnp.ndarray, usm_ndarray} + Matrices for which the eigenvalues and right eigenvectors will + be computed. + + Returns + ------- + eigenvalues : (..., M) dpnp.ndarray + The eigenvalues, each repeated according to its multiplicity. + The eigenvalues are not necessarily ordered. The resulting + array will be of complex type, unless the imaginary part is + zero in which case it will be cast to a real type. When `a` + is real the resulting eigenvalues will be real (0 imaginary + part) or occur in conjugate pairs + eigenvectors : (..., M, M) dpnp.ndarray + The normalized (unit "length") eigenvectors, such that the + column ``v[:,i]`` is the eigenvector corresponding to the + eigenvalue ``w[i]``. + + Note + ---- + Since there is no proper OneMKL LAPACK function, DPNP will calculate + through a fallback on NumPy call. + + See Also + -------- + :obj:`dpnp.linalg.eigvals` : Compute the eigenvalues of a general matrix. + :obj:`dpnp.linalg.eigh` : Return the eigenvalues and eigenvectors of a complex Hermitian + (conjugate symmetric) or a real symmetric matrix. + :obj:`dpnp.linalg.eigvalsh` : Compute the eigenvalues of a complex Hermitian or + real symmetric matrix. + + Examples + -------- + >>> import dpnp as np + >>> from dpnp import linalg as LA + + (Almost) trivial example with real eigenvalues and eigenvectors. + + >>> w, v = LA.eig(np.diag((1, 2, 3))) + >>> w, v + (array([1., 2., 3.]), + array([[1., 0., 0.], + [0., 1., 0.], + [0., 0., 1.]])) + + Real matrix possessing complex eigenvalues and eigenvectors; + note that the eigenvalues are complex conjugates of each other. + + >>> w, v = LA.eig(np.array([[1, -1], [1, 1]])) + >>> w, v + (array([1.+1.j, 1.-1.j]), + array([[0.70710678+0.j , 0.70710678-0.j ], + [0. -0.70710678j, 0. +0.70710678j]])) + + Complex-valued matrix with real eigenvalues (but complex-valued + eigenvectors); note that ``a.conj().T == a``, i.e., `a` is Hermitian. + + >>> a = np.array([[1, 1j], [-1j, 1]]) + >>> w, v = LA.eig(a) + >>> w, v + (array([2.+0.j, 0.+0.j]), + array([[ 0. +0.70710678j, 0.70710678+0.j ], # may vary + [ 0.70710678+0.j , -0. +0.70710678j]]) + + Be careful about round-off error! + + >>> a = np.array([[1 + 1e-9, 0], [0, 1 - 1e-9]]) + >>> # Theor. eigenvalues are 1 +/- 1e-9 + >>> w, v = LA.eig(a) + >>> w, v + (array([1., 1.]), + array([[1., 0.], + [0., 1.]])) """ - x1_desc = dpnp.get_dpnp_descriptor(x1, copy_when_nondefault_queue=False) - if x1_desc: - if x1_desc.size > 0: - return dpnp_eig(x1_desc) + dpnp.check_supported_arrays_type(a) + check_stacked_2d(a) + check_stacked_square(a) - return call_origin(numpy.linalg.eig, x1) + a_sycl_queue = a.sycl_queue + a_usm_type = a.usm_type + + # Since geev function from OneMKL LAPACK is not implemented yet, + # use NumPy for this calculation. + w_np, v_np = numpy.linalg.eig(dpnp.asnumpy(a)) + return ( + dpnp.array(w_np, sycl_queue=a_sycl_queue, usm_type=a_usm_type), + dpnp.array(v_np, sycl_queue=a_sycl_queue, usm_type=a_usm_type), + ) def eigh(a, UPLO="L"): @@ -334,32 +415,73 @@ def eigh(a, UPLO="L"): return dpnp_eigh(a, UPLO=UPLO) -def eigvals(input): +def eigvals(a): """ Compute the eigenvalues of a general matrix. - Main difference between `eigvals` and `eig`: the eigenvectors aren't - returned. + For full documentation refer to :obj:`numpy.linalg.eigvals`. Parameters ---------- - input : (..., M, M) array_like + a : (..., M, M) {dpnp.ndarray, usm_ndarray} A complex- or real-valued matrix whose eigenvalues will be computed. Returns ------- - w : (..., M,) ndarray + w : (..., M) dpnp.ndarray The eigenvalues, each repeated according to its multiplicity. They are not necessarily ordered, nor are they necessarily real for real matrices. + + Note + ---- + Since there is no proper OneMKL LAPACK function, DPNP will calculate + through a fallback on NumPy call. + + See Also + -------- + :obj:`dpnp.linalg.eig` : Compute the eigenvalues and right eigenvectors of a square array. + :obj:`dpnp.linalg.eigvalsh` : Compute the eigenvalues of a complex Hermitian or + real symmetric matrix. + :obj:`dpnp.linalg.eigh` : Return the eigenvalues and eigenvectors of a complex Hermitian + (conjugate symmetric) or a real symmetric matrix. + + Examples + -------- + Illustration, using the fact that the eigenvalues of a diagonal matrix + are its diagonal elements, that multiplying a matrix on the left + by an orthogonal matrix, `Q`, and on the right by `Q.T` (the transpose + of `Q`), preserves the eigenvalues of the "middle" matrix. In other words, + if `Q` is orthogonal, then ``Q * A * Q.T`` has the same eigenvalues as + ``A``: + + >>> import dpnp as np + >>> from dpnp import linalg as LA + >>> x = np.random.random() + >>> Q = np.array([[np.cos(x), -np.sin(x)], [np.sin(x), np.cos(x)]]) + >>> LA.norm(Q[0, :]), LA.norm(Q[1, :]), np.dot(Q[0, :],Q[1, :]) + (array(1.), array(1.), array(0.)) + + Now multiply a diagonal matrix by ``Q`` on one side and by ``Q.T`` on the other: + + >>> D = np.diag((-1,1)) + >>> LA.eigvals(D) + array([-1., 1.]) + >>> A = np.dot(Q, D) + >>> A = np.dot(A, Q.T) + >>> LA.eigvals(A) + array([-1., 1.]) # random + """ - x1_desc = dpnp.get_dpnp_descriptor(input, copy_when_nondefault_queue=False) - if x1_desc: - if x1_desc.size > 0: - return dpnp_eigvals(x1_desc).get_pyobj() + dpnp.check_supported_arrays_type(a) + check_stacked_2d(a) + check_stacked_square(a) - return call_origin(numpy.linalg.eigvals, input) + # Since geev function from OneMKL LAPACK is not implemented yet, + # use NumPy for this calculation. + w_np = numpy.linalg.eigvals(dpnp.asnumpy(a)) + return dpnp.array(w_np, sycl_queue=a.sycl_queue, usm_type=a.usm_type) def eigvalsh(a, UPLO="L"): @@ -493,6 +615,8 @@ def matrix_power(a, n): elements is the same as those of `M`. If the exponent is negative the elements are floating-point. + Examples + -------- >>> import dpnp as np >>> i = np.array([[0, 1], [-1, 0]]) # matrix equiv. of the imaginary unit >>> np.linalg.matrix_power(i, 3) # should = -i @@ -881,16 +1005,17 @@ def solve(a, b): For full documentation refer to :obj:`numpy.linalg.solve`. + Parameters + ---------- + a : (..., M, M) {dpnp.ndarray, usm_ndarray} + Coefficient matrix. + b : {(…, M,), (…, M, K)} {dpnp.ndarray, usm_ndarray} + Ordinate or "dependent variable" values. + Returns ------- out : {(…, M,), (…, M, K)} dpnp.ndarray - Solution to the system ax = b. Returned shape is identical to b. - - Limitations - ----------- - 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`. + Solution to the system `ax = b`. Returned shape is identical to `b`. See Also -------- diff --git a/tests/test_linalg.py b/tests/test_linalg.py index e2ade89d9ea..92f70c70215 100644 --- a/tests/test_linalg.py +++ b/tests/test_linalg.py @@ -456,48 +456,6 @@ def test_det_errors(self): assert_raises(TypeError, inp.linalg.det, a_np) -@pytest.mark.parametrize("type", get_all_dtypes(no_bool=True, no_complex=True)) -@pytest.mark.parametrize("size", [2, 4, 8, 16, 300]) -def test_eig_arange(type, size): - a = numpy.arange(size * size, dtype=type).reshape((size, size)) - symm_orig = ( - numpy.tril(a) - + numpy.tril(a, -1).T - + numpy.diag(numpy.full((size,), size * size, dtype=type)) - ) - symm = symm_orig - dpnp_symm_orig = inp.array(symm) - dpnp_symm = dpnp_symm_orig - - dpnp_val, dpnp_vec = inp.linalg.eig(dpnp_symm) - np_val, np_vec = numpy.linalg.eig(symm) - - # DPNP sort val/vec by abs value - vvsort(dpnp_val, dpnp_vec, size, inp) - - # NP sort val/vec by abs value - vvsort(np_val, np_vec, size, numpy) - - # NP change sign of vectors - for i in range(np_vec.shape[1]): - if np_vec[0, i] * dpnp_vec[0, i] < 0: - np_vec[:, i] = -np_vec[:, i] - - assert_array_equal(symm_orig, symm) - assert_array_equal(dpnp_symm_orig, dpnp_symm) - - if has_support_aspect64(): - assert dpnp_val.dtype == np_val.dtype - assert dpnp_vec.dtype == np_vec.dtype - assert dpnp_val.shape == np_val.shape - assert dpnp_vec.shape == np_vec.shape - assert dpnp_val.usm_type == dpnp_symm.usm_type - assert dpnp_vec.usm_type == dpnp_symm.usm_type - - assert_allclose(dpnp_val, np_val, rtol=1e-05, atol=1e-05) - assert_allclose(dpnp_vec, np_vec, rtol=1e-05, atol=1e-05) - - class TestEigenvalue: # Eigenvalue decomposition of a matrix or a batch of matrices # by checking if the eigen equation A*v=w*v holds for given eigenvalues(w) @@ -519,6 +477,8 @@ def assert_eigen_decomposition(self, a, w, v, rtol=1e-5, atol=1e-5): @pytest.mark.parametrize( "func", [ + "eig", + "eigvals", "eigh", "eigvalsh", ], @@ -537,8 +497,12 @@ def assert_eigen_decomposition(self, a, w, v, rtol=1e-5, atol=1e-5): ], ) def test_eigenvalues(self, func, shape, dtype, order): + # Set a `hermitian` flag for generate_random_numpy_array() to + # get a symmetric array for eigh() and eigvalsh() or + # non-symmetric for eig() and eigvals() + is_hermitian = func in ("eigh, eigvalsh") a = generate_random_numpy_array( - shape, dtype, hermitian=True, seed_value=81 + shape, dtype, hermitian=is_hermitian, seed_value=81 ) a_order = numpy.array(a, order=order) a_dp = inp.array(a, order=order) @@ -548,21 +512,53 @@ def test_eigenvalues(self, func, shape, dtype, order): # However, both OneMKL and rocSOLVER pick a different convention for # constructing eigenvectors, so v's are not directly comparible and # we verify them through the eigen equation A*v=w*v. - if func == "eigh": - w, _ = numpy.linalg.eigh(a_order) - w_dp, v_dp = inp.linalg.eigh(a_dp) + if func in ("eig", "eigh"): + w, _ = getattr(numpy.linalg, func)(a_order) + w_dp, v_dp = getattr(inp.linalg, func)(a_dp) self.assert_eigen_decomposition(a_dp, w_dp, v_dp) - else: # eighvalsh - w = numpy.linalg.eigvalsh(a_order) - w_dp = inp.linalg.eigvalsh(a_dp) + else: # eighvals or eigvalsh + w = getattr(numpy.linalg, func)(a_order) + w_dp = getattr(inp.linalg, func)(a_dp) assert_dtype_allclose(w_dp, w) + # eigh() and eigvalsh() are tested in cupy tests @pytest.mark.parametrize( "func", [ + "eig", + "eigvals", + ], + ) + @pytest.mark.parametrize( + "shape", + [(0, 0), (2, 0, 0), (0, 3, 3)], + ids=["(0,0)", "(2,0,0)", "(0,3,3)"], + ) + @pytest.mark.parametrize("dtype", get_all_dtypes(no_bool=True)) + def test_eigenvalue_empty(self, func, shape, dtype): + a_np = numpy.empty(shape, dtype=dtype) + a_dp = inp.array(a_np) + + if func == "eig": + w, v = getattr(numpy.linalg, func)(a_np) + w_dp, v_dp = getattr(inp.linalg, func)(a_dp) + + assert_dtype_allclose(v_dp, v) + + else: # eigvals + w = getattr(numpy.linalg, func)(a_np) + w_dp = getattr(inp.linalg, func)(a_dp) + + assert_dtype_allclose(w_dp, w) + + @pytest.mark.parametrize( + "func", + [ + "eig", + "eigvals", "eigh", "eigvalsh", ], @@ -584,22 +580,8 @@ def test_eigenvalue_errors(self, func): assert_raises(inp.linalg.LinAlgError, dpnp_func, a_dp_not_scquare) # invalid UPLO - assert_raises(ValueError, dpnp_func, a_dp, UPLO="N") - - -@pytest.mark.parametrize("type", get_all_dtypes(no_bool=True, no_complex=True)) -def test_eigvals(type): - if dpctl.SyclDevice().device_type != dpctl.device_type.gpu: - pytest.skip( - "eigvals function doesn't work on CPU: https://github.com/IntelPython/dpnp/issues/1005" - ) - arrays = [[[0, 0], [0, 0]], [[1, 2], [1, 2]], [[1, 2], [3, 4]]] - for array in arrays: - a = numpy.array(array, dtype=type) - ia = inp.array(a) - result = inp.linalg.eigvals(ia) - expected = numpy.linalg.eigvals(a) - assert_allclose(expected, result, atol=0.5) + if func in ("eigh", "eigvalsh"): + assert_raises(ValueError, dpnp_func, a_dp, UPLO="N") class TestInv: diff --git a/tests/test_sycl_queue.py b/tests/test_sycl_queue.py index 2b14bbb330f..08db52ad70c 100644 --- a/tests/test_sycl_queue.py +++ b/tests/test_sycl_queue.py @@ -1189,64 +1189,11 @@ def test_det(device): assert_sycl_queue_equal(result_queue, expected_queue) -@pytest.mark.usefixtures("allow_fall_back_on_numpy") -@pytest.mark.parametrize( - "device", - valid_devices, - ids=[device.filter_string for device in valid_devices], -) -def test_eig(device): - if device.device_type != dpctl.device_type.gpu: - pytest.skip( - "eig function doesn't work on CPU: https://github.com/IntelPython/dpnp/issues/1005" - ) - - size = 4 - dtype = dpnp.default_float_type(device) - a = numpy.arange(size * size, dtype=dtype).reshape((size, size)) - symm_orig = ( - numpy.tril(a) - + numpy.tril(a, -1).T - + numpy.diag(numpy.full((size,), size * size, dtype=dtype)) - ) - numpy_data = symm_orig - dpnp_symm_orig = dpnp.array(numpy_data, device=device) - dpnp_data = dpnp_symm_orig - - dpnp_val, dpnp_vec = dpnp.linalg.eig(dpnp_data) - numpy_val, numpy_vec = numpy.linalg.eig(numpy_data) - - # DPNP sort val/vec by abs value - vvsort(dpnp_val, dpnp_vec, size, dpnp) - - # NP sort val/vec by abs value - vvsort(numpy_val, numpy_vec, size, numpy) - - # NP change sign of vectors - for i in range(numpy_vec.shape[1]): - if numpy_vec[0, i] * dpnp_vec[0, i] < 0: - numpy_vec[:, i] = -numpy_vec[:, i] - - assert_allclose(dpnp_val, numpy_val, rtol=1e-05, atol=1e-05) - assert_allclose(dpnp_vec, numpy_vec, rtol=1e-05, atol=1e-05) - - assert dpnp_val.dtype == numpy_val.dtype - assert dpnp_vec.dtype == numpy_vec.dtype - assert dpnp_val.shape == numpy_val.shape - assert dpnp_vec.shape == numpy_vec.shape - - expected_queue = dpnp_data.get_array().sycl_queue - dpnp_val_queue = dpnp_val.get_array().sycl_queue - dpnp_vec_queue = dpnp_vec.get_array().sycl_queue - - # compare queue and device - assert_sycl_queue_equal(dpnp_val_queue, expected_queue) - assert_sycl_queue_equal(dpnp_vec_queue, expected_queue) - - @pytest.mark.parametrize( "func", [ + "eig", + "eigvals", "eigh", "eigvalsh", ], @@ -1275,16 +1222,22 @@ def test_eig(device): ) def test_eigenvalue(func, shape, device): dtype = dpnp.default_float_type(device) + # Set a `hermitian` flag for generate_random_numpy_array() to + # get a symmetric array for eigh() and eigvalsh() or + # non-symmetric for eig() and eigvals() + is_hermitian = func in ("eigh, eigvalsh") # Set seed_value=81 to prevent # random generation of the input singular matrix - a = generate_random_numpy_array(shape, dtype, hermitian=True, seed_value=81) + a = generate_random_numpy_array( + shape, dtype, hermitian=is_hermitian, seed_value=81 + ) dp_a = dpnp.array(a, device=device) expected_queue = dp_a.get_array().sycl_queue - if func == "eigh": - dp_val, dp_vec = dpnp.linalg.eigh(dp_a) - np_val, np_vec = numpy.linalg.eigh(a) + if func in ("eig", "eigh"): + dp_val, dp_vec = getattr(dpnp.linalg, func)(dp_a) + np_val, np_vec = getattr(numpy.linalg, func)(a) # Check the eigenvalue decomposition if a.ndim == 2: @@ -1306,9 +1259,9 @@ def test_eigenvalue(func, shape, device): # compare queue and device assert_sycl_queue_equal(dpnp_vec_queue, expected_queue) - else: # eighvalsh - dp_val = dpnp.linalg.eigvalsh(dp_a) - np_val = numpy.linalg.eigvalsh(a) + else: # eighvals or eigvalsh + dp_val = getattr(dpnp.linalg, func)(dp_a) + np_val = getattr(numpy.linalg, func)(a) assert_allclose(dp_val, np_val, rtol=1e-05, atol=1e-05) assert dp_val.shape == np_val.shape @@ -1319,31 +1272,6 @@ def test_eigenvalue(func, shape, device): assert_sycl_queue_equal(dpnp_val_queue, expected_queue) -@pytest.mark.parametrize( - "device", - valid_devices, - ids=[device.filter_string for device in valid_devices], -) -def test_eigvals(device): - if device.device_type != dpctl.device_type.gpu: - pytest.skip( - "eigvals function doesn't work on CPU: https://github.com/IntelPython/dpnp/issues/1005" - ) - - data = [[0, 0], [0, 0]] - numpy_data = numpy.array(data) - dpnp_data = dpnp.array(data, device=device) - - result = dpnp.linalg.eigvals(dpnp_data) - expected = numpy.linalg.eigvals(numpy_data) - assert_allclose(expected, result, atol=0.5) - - expected_queue = dpnp_data.get_array().sycl_queue - result_queue = result.get_array().sycl_queue - - assert_sycl_queue_equal(result_queue, expected_queue) - - @pytest.mark.parametrize( "shape, is_empty", [ diff --git a/tests/test_usm_type.py b/tests/test_usm_type.py index 3b2aad98ea2..f95c0413053 100644 --- a/tests/test_usm_type.py +++ b/tests/test_usm_type.py @@ -829,6 +829,8 @@ def test_where(usm_type): @pytest.mark.parametrize( "func", [ + "eig", + "eigvals", "eigh", "eigvalsh", ], @@ -852,15 +854,19 @@ def test_where(usm_type): ) @pytest.mark.parametrize("usm_type", list_of_usm_types, ids=list_of_usm_types) def test_eigenvalue(func, shape, usm_type): - a_np = generate_random_numpy_array(shape, hermitian=True) + # Set a `hermitian` flag for generate_random_numpy_array() to + # get a symmetric array for eigh() and eigvalsh() or + # non-symmetric for eig() and eigvals() + is_hermitian = func in ("eigh, eigvalsh") + a_np = generate_random_numpy_array(shape, hermitian=is_hermitian) a = dp.array(a_np, usm_type=usm_type) - if func == "eigh": - dp_val, dp_vec = dp.linalg.eigh(a) + if func in ("eig", "eigh"): + dp_val, dp_vec = getattr(dp.linalg, func)(a) assert a.usm_type == dp_vec.usm_type - else: # eighvalsh - dp_val = dp.linalg.eigvalsh(a) + else: # eighvals or eigvalsh + dp_val = getattr(dp.linalg, func)(a) assert a.usm_type == dp_val.usm_type