From 4c283948d456e5affe4cb76452cffbd1217c341e Mon Sep 17 00:00:00 2001 From: Vladislav Perevezentsev Date: Thu, 15 Feb 2024 14:11:22 +0100 Subject: [PATCH 01/26] Add dpnp.linalg.eigvalsh() impl --- dpnp/linalg/dpnp_iface_linalg.py | 59 ++++++++++++++++++++++++++++++++ dpnp/linalg/dpnp_utils_linalg.py | 27 ++++++++++++--- 2 files changed, 82 insertions(+), 4 deletions(-) diff --git a/dpnp/linalg/dpnp_iface_linalg.py b/dpnp/linalg/dpnp_iface_linalg.py index 88a904b3c3c4..5b736ee09cda 100644 --- a/dpnp/linalg/dpnp_iface_linalg.py +++ b/dpnp/linalg/dpnp_iface_linalg.py @@ -50,6 +50,7 @@ dpnp_cholesky, dpnp_det, dpnp_eigh, + dpnp_eigvalsh, dpnp_inv, dpnp_qr, dpnp_slogdet, @@ -64,6 +65,7 @@ "eig", "eigh", "eigvals", + "eigvalsh", "inv", "matrix_power", "matrix_rank", @@ -319,6 +321,63 @@ def eigvals(input): return call_origin(numpy.linalg.eigvals, input) +def eigvalsh(a, UPLO="L"): + """ + eigvalsh(a, UPLO="L") + + Return the eigenvalues of a complex Hermitian or real symmetric matrix. + + Main difference from :obj:`dpnp.linalg.eigh`: the eigenvectors are not computed. + + For full documentation refer to :obj:`numpy.linalg.eigvalsh`. + + Parameters + ---------- + a : (..., M) {dpnp.ndarray, usm_ndarray} + A complex- or real-valued array whose eigenvalues are to be computed. + UPLO : {"L", "U"}, optional + Specifies the calculation uses either the lower ("L") or upper ("U") + triangular part of the matrix. + Regardless of this choice, only the real parts of the diagonal are + considered to preserve the Hermite matrix property. + It therefore follows that the imaginary part of the diagonal + will always be treated as zero. + Default: "L". + + Returns + ------- + w : (..., M) dpnp.ndarray + The eigenvalues in ascending order, each repeated according to + its multiplicity. + + See Also + -------- + :obj:`dpnp.linalg.eigh` : Returns the eigenvalues and eigenvectors of real symmetric + or complex Hermitian (conjugate symmetric) arrays. + :obj:`dpnp.linalg.eigvals` : Returns the eigenvalues of general real or complex arrays. + :obj:`dpnp.linalg.eig` : Returns the eigenvalues and right eigenvectors of + general real or complex arrays. + + Examples + -------- + >>> import dpnp as np + >>> from dpnp import linalg as LA + >>> a = np.array([[1, -2j], [2j, 5]]) + >>> LA.eigvalsh(a) + array([0.17157288, 5.82842712]) + + """ + + dpnp.check_supported_arrays_type(a) + check_stacked_2d(a) + check_stacked_square(a) + + if UPLO not in ("L", "U"): + raise ValueError("UPLO argument must be 'L' or 'U'") + + return dpnp_eigvalsh(a, UPLO=UPLO) + + def inv(a): """ Compute the (multiplicative) inverse of a matrix. diff --git a/dpnp/linalg/dpnp_utils_linalg.py b/dpnp/linalg/dpnp_utils_linalg.py index a6dcfbf0c2bf..7e1dd8084f62 100644 --- a/dpnp/linalg/dpnp_utils_linalg.py +++ b/dpnp/linalg/dpnp_utils_linalg.py @@ -38,6 +38,7 @@ "dpnp_cholesky", "dpnp_det", "dpnp_eigh", + "dpnp_eigvalsh", "dpnp_inv", "dpnp_qr", "dpnp_slogdet", @@ -731,9 +732,9 @@ def dpnp_det(a): return det.reshape(shape) -def dpnp_eigh(a, UPLO): +def dpnp_eigh(a, UPLO, eigen_mode="V"): """ - dpnp_eigh(a, UPLO) + dpnp_eigh(a, UPLO, eigen_mode="V") Return the eigenvalues and eigenvectors of a complex Hermitian (conjugate symmetric) or a real symmetric matrix. @@ -749,8 +750,11 @@ def dpnp_eigh(a, UPLO): a_order = "C" if a.flags.c_contiguous else "F" a_usm_arr = dpnp.get_usm_ndarray(a) - # 'V' means both eigenvectors and eigenvalues will be calculated - jobz = _jobz["V"] + # `eigen_mode` can be either "N" or "V", specifying the computation mode + # for OneMKL LAPACK `syevd` and `heevd` routines. + # "V" (default) means both eigenvectors and eigenvalues will be calculated + # "N" means only eigenvalues will be calculated + jobz = _jobz[eigen_mode] uplo = _upper_lower[UPLO] # get resulting type of arrays with eigenvalues and eigenvectors @@ -859,6 +863,21 @@ def dpnp_eigh(a, UPLO): return w, out_v +def dpnp_eigvalsh(a, UPLO): + """ + dpnp_eigvalsh(a, UPLO) + + Return the eigenvalues of a complex Hermitian or real symmetric matrix. + + """ + + if a.size == 0: + res_type = _real_type(_common_type(a)) + return dpnp.empty_like(a, shape=a.shape[:-1], dtype=res_type) + + return dpnp_eigh(a, UPLO=UPLO, eigen_mode="N")[0] + + def dpnp_inv_batched(a, res_type): """ dpnp_inv_batched(a, res_type) From 0c891efa100396fad7a1a9acc7e99f1f059c46d2 Mon Sep 17 00:00:00 2001 From: Vladislav Perevezentsev Date: Thu, 15 Feb 2024 14:12:33 +0100 Subject: [PATCH 02/26] Update cupy tests for dpnp.linalg.eigvalsh func --- .../cupy/linalg_tests/test_eigenvalue.py | 30 ++++++++++--------- 1 file changed, 16 insertions(+), 14 deletions(-) diff --git a/tests/third_party/cupy/linalg_tests/test_eigenvalue.py b/tests/third_party/cupy/linalg_tests/test_eigenvalue.py index b620bd39e984..93e860e5a9f1 100644 --- a/tests/third_party/cupy/linalg_tests/test_eigenvalue.py +++ b/tests/third_party/cupy/linalg_tests/test_eigenvalue.py @@ -113,9 +113,10 @@ def test_eigh_complex_batched(self, xp, dtype): ) return w - @pytest.mark.skip("No support of dpnp.eigvalsh()") @testing.for_all_dtypes(no_float16=True, no_complex=True) - @testing.numpy_cupy_allclose(rtol=1e-3, atol=1e-4) + @testing.numpy_cupy_allclose( + rtol=1e-3, atol=1e-4, type_check=has_support_aspect64() + ) def test_eigvalsh(self, xp, dtype): a = xp.array([[1, 0, 3], [0, 5, 0], [7, 0, 9]], dtype) w = xp.linalg.eigvalsh(a, UPLO=self.UPLO) @@ -123,9 +124,10 @@ def test_eigvalsh(self, xp, dtype): # so they should be directly comparable return w - @pytest.mark.skip("No support of dpnp.eigvalsh()") @testing.for_all_dtypes(no_float16=True, no_complex=True) - @testing.numpy_cupy_allclose(rtol=1e-3, atol=1e-4) + @testing.numpy_cupy_allclose( + rtol=1e-3, atol=1e-4, type_check=has_support_aspect64() + ) def test_eigvalsh_batched(self, xp, dtype): a = xp.array( [ @@ -139,9 +141,10 @@ def test_eigvalsh_batched(self, xp, dtype): # so they should be directly comparable return w - @pytest.mark.skip("No support of dpnp.eigvalsh()") @testing.for_complex_dtypes() - @testing.numpy_cupy_allclose(rtol=1e-3, atol=1e-4) + @testing.numpy_cupy_allclose( + rtol=1e-3, atol=1e-4, type_check=has_support_aspect64() + ) def test_eigvalsh_complex(self, xp, dtype): a = xp.array([[1, 2j, 3], [4j, 5, 6j], [7, 8j, 9]], dtype) w = xp.linalg.eigvalsh(a, UPLO=self.UPLO) @@ -149,9 +152,10 @@ def test_eigvalsh_complex(self, xp, dtype): # so they should be directly comparable return w - @pytest.mark.skip("No support of dpnp.eigvalsh()") @testing.for_complex_dtypes() - @testing.numpy_cupy_allclose(rtol=1e-3, atol=1e-4) + @testing.numpy_cupy_allclose( + rtol=1e-3, atol=1e-4, type_check=has_support_aspect64() + ) def test_eigvalsh_complex_batched(self, xp, dtype): a = xp.array( [ @@ -179,11 +183,10 @@ def test_eigh(self, xp, dtype): assert a.size == 0 return xp.linalg.eigh(a, UPLO=self.UPLO) - @pytest.mark.skip("No support of dpnp.eigvalsh()") @testing.for_dtypes("ifdFD") - @testing.numpy_cupy_allclose() + @testing.numpy_cupy_allclose(type_check=has_support_aspect64()) def test_eigvalsh(self, xp, dtype): - a = xp.empty(self.shape, dtype) + a = xp.empty(self.shape, dtype=dtype) assert a.size == 0 return xp.linalg.eigvalsh(a, UPLO=self.UPLO) @@ -192,7 +195,7 @@ def test_eigvalsh(self, xp, dtype): *testing.product( { "UPLO": ["U", "L"], - "shape": [(), (3,), (2, 3), (4, 0), (2, 2, 3), (0, 2, 3)], + "shape": [()], } ) ) @@ -203,9 +206,8 @@ def test_eigh_shape_error(self): with pytest.raises((numpy.linalg.LinAlgError, ValueError)): xp.linalg.eigh(a, self.UPLO) - @pytest.mark.skip("No support of dpnp.eigvalsh()") def test_eigvalsh_shape_error(self): for xp in (numpy, cupy): a = xp.zeros(self.shape) - with pytest.raises((numpy.linalg.LinAlgError, ValueError)): + with pytest.raises(xp.linalg.LinAlgError): xp.linalg.eigvalsh(a, self.UPLO) From 5e65e477579e0d4ab566c10a8c8d0449f543f54f Mon Sep 17 00:00:00 2001 From: Vladislav Perevezentsev Date: Thu, 15 Feb 2024 15:22:16 +0100 Subject: [PATCH 03/26] Small update dpnp.linalg.eigh --- dpnp/linalg/dpnp_iface_linalg.py | 47 +++++++++++++-------------- dpnp/linalg/dpnp_utils_linalg.py | 54 +++++++++++++------------------- 2 files changed, 46 insertions(+), 55 deletions(-) diff --git a/dpnp/linalg/dpnp_iface_linalg.py b/dpnp/linalg/dpnp_iface_linalg.py index 5b736ee09cda..945d03df4c6f 100644 --- a/dpnp/linalg/dpnp_iface_linalg.py +++ b/dpnp/linalg/dpnp_iface_linalg.py @@ -241,6 +241,19 @@ def eigh(a, UPLO="L"): For full documentation refer to :obj:`numpy.linalg.eigh`. + Parameters + ---------- + a : (..., M, M) {dpnp.ndarray, usm_ndarray} + A complex- or real-valued array whose eigenvalues and eigenvectors are to be computed. + UPLO : {"L", "U"}, optional + Specifies the calculation uses either the lower ("L") or upper ("U") + triangular part of the matrix. + Regardless of this choice, only the real parts of the diagonal are + considered to preserve the Hermite matrix property. + It therefore follows that the imaginary part of the diagonal + will always be treated as zero. + Default: "L". + Returns ------- w : (..., M) dpnp.ndarray @@ -250,15 +263,11 @@ def eigh(a, UPLO="L"): The column ``v[:, i]`` is the normalized eigenvector corresponding to the eigenvalue ``w[i]``. - Limitations - ----------- - Parameter `a` is supported as :class:`dpnp.ndarray` or :class:`dpctl.tensor.usm_ndarray`. - Input array data types are limited by supported DPNP :ref:`Data types`. - See Also -------- - :obj:`dpnp.eig` : eigenvalues and right eigenvectors for non-symmetric arrays. - :obj:`dpnp.eigvals` : eigenvalues of non-symmetric arrays. + :obj:`dpnp.linalg.eig` : Return the eigenvalues and right eigenvectors of + a general matrix. + :obj:`dpnp.linalg.eigvals` : Return the eigenvalues of a general matrix. Examples -------- @@ -276,20 +285,12 @@ def eigh(a, UPLO="L"): """ dpnp.check_supported_arrays_type(a) + check_stacked_2d(a) + check_stacked_square(a) if UPLO not in ("L", "U"): raise ValueError("UPLO argument must be 'L' or 'U'") - if a.ndim < 2: - raise ValueError( - "%d-dimensional array given. Array must be " - "at least two-dimensional" % a.ndim - ) - - m, n = a.shape[-2:] - if m != n: - raise ValueError("Last 2 dimensions of the array must be square") - return dpnp_eigh(a, UPLO=UPLO) @@ -333,7 +334,7 @@ def eigvalsh(a, UPLO="L"): Parameters ---------- - a : (..., M) {dpnp.ndarray, usm_ndarray} + a : (..., M, M) {dpnp.ndarray, usm_ndarray} A complex- or real-valued array whose eigenvalues are to be computed. UPLO : {"L", "U"}, optional Specifies the calculation uses either the lower ("L") or upper ("U") @@ -352,11 +353,11 @@ def eigvalsh(a, UPLO="L"): See Also -------- - :obj:`dpnp.linalg.eigh` : Returns the eigenvalues and eigenvectors of real symmetric - or complex Hermitian (conjugate symmetric) arrays. - :obj:`dpnp.linalg.eigvals` : Returns the eigenvalues of general real or complex arrays. - :obj:`dpnp.linalg.eig` : Returns the eigenvalues and right eigenvectors of - general real or complex arrays. + :obj:`dpnp.linalg.eigh` : Return the eigenvalues and eigenvectors of a complex Hermitian + (conjugate symmetric) or a real symmetric matrix. + :obj:`dpnp.linalg.eigvals` : Return the eigenvalues of a general matrix. + :obj:`dpnp.linalg.eig` : Return the eigenvalues and right eigenvectors of + a general matrix. Examples -------- diff --git a/dpnp/linalg/dpnp_utils_linalg.py b/dpnp/linalg/dpnp_utils_linalg.py index 7e1dd8084f62..676d0bc886ab 100644 --- a/dpnp/linalg/dpnp_utils_linalg.py +++ b/dpnp/linalg/dpnp_utils_linalg.py @@ -745,7 +745,6 @@ def dpnp_eigh(a, UPLO, eigen_mode="V"): """ - a_usm_type = a.usm_type a_sycl_queue = a.sycl_queue a_order = "C" if a.flags.c_contiguous else "F" a_usm_arr = dpnp.get_usm_ndarray(a) @@ -758,49 +757,43 @@ def dpnp_eigh(a, UPLO, eigen_mode="V"): uplo = _upper_lower[UPLO] # get resulting type of arrays with eigenvalues and eigenvectors - a_dtype = a.dtype - lapack_func = "_syevd" - if dpnp.issubdtype(a_dtype, dpnp.complexfloating): - lapack_func = "_heevd" - v_type = a_dtype - w_type = dpnp.float64 if a_dtype == dpnp.complex128 else dpnp.float32 - elif dpnp.issubdtype(a_dtype, dpnp.floating): - v_type = w_type = a_dtype - elif a_sycl_queue.sycl_device.has_aspect_fp64: - v_type = w_type = dpnp.float64 - else: - v_type = w_type = dpnp.float32 + v_type = _common_type(a) + w_type = _real_type(v_type) + + # Get LAPACK function (_syevd for real or _heevd for complex data types) + # to compute all eigenvalues and, optionally, all eigenvectors + lapack_func = ( + "_heevd" if dpnp.issubdtype(v_type, dpnp.complexfloating) else "_syevd" + ) if a.ndim > 2: - w = dpnp.empty( - a.shape[:-1], + w = dpnp.empty_like( + a, + shape=a.shape[:-1], dtype=w_type, - usm_type=a_usm_type, - sycl_queue=a_sycl_queue, ) # need to loop over the 1st dimension to get eigenvalues and eigenvectors of 3d matrix A - op_count = a.shape[0] - if op_count == 0: + batch_size = a.shape[0] + if batch_size == 0: return w, dpnp.empty_like(a, dtype=v_type) - eig_vecs = [None] * op_count - ht_copy_ev = [None] * op_count - ht_lapack_ev = [None] * op_count - for i in range(op_count): + eig_vecs = [None] * batch_size + ht_list_ev = [None] * batch_size * 2 + 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 eigenvectors eig_vecs[i] = dpnp.empty_like(a[i], order="F", dtype=v_type) # use DPCTL tensor function to fill the array of eigenvectors with content of input array - ht_copy_ev[i], copy_ev = ti._copy_usm_ndarray_into_usm_ndarray( + ht_list_ev[2 * i], copy_ev = ti._copy_usm_ndarray_into_usm_ndarray( src=a_usm_arr[i], dst=eig_vecs[i].get_array(), sycl_queue=a_sycl_queue, ) # call LAPACK extension function to get eigenvalues and eigenvectors of a portion of matrix A - ht_lapack_ev[i], _ = getattr(li, lapack_func)( + ht_list_ev[2 * i + 1], _ = getattr(li, lapack_func)( a_sycl_queue, jobz, uplo, @@ -809,9 +802,7 @@ def dpnp_eigh(a, UPLO, eigen_mode="V"): depends=[copy_ev], ) - for i in range(op_count): - ht_lapack_ev[i].wait() - ht_copy_ev[i].wait() + dpctl.SyclEvent.wait_for(ht_list_ev) # combine the list of eigenvectors into a single array v = dpnp.array(eig_vecs, order=a_order) @@ -827,11 +818,10 @@ def dpnp_eigh(a, UPLO, eigen_mode="V"): ) # allocate a memory for dpnp array of eigenvalues - w = dpnp.empty( - a.shape[:-1], + w = dpnp.empty_like( + a, + shape=a.shape[:-1], dtype=w_type, - usm_type=a_usm_type, - sycl_queue=a_sycl_queue, ) # call LAPACK extension function to get eigenvalues and eigenvectors of matrix A From 5c6a7ceef2842c29de88e4162e0f3c6d2e464d59 Mon Sep 17 00:00:00 2001 From: Vladislav Perevezentsev Date: Thu, 15 Feb 2024 15:29:29 +0100 Subject: [PATCH 04/26] Update cupy tests for dpnp.linalg.eigh --- tests/skipped_tests_gpu_no_fp64.tbl | 3 --- .../cupy/linalg_tests/test_eigenvalue.py | 13 ++++++++----- 2 files changed, 8 insertions(+), 8 deletions(-) diff --git a/tests/skipped_tests_gpu_no_fp64.tbl b/tests/skipped_tests_gpu_no_fp64.tbl index d724a6043e57..56045d9eb52b 100644 --- a/tests/skipped_tests_gpu_no_fp64.tbl +++ b/tests/skipped_tests_gpu_no_fp64.tbl @@ -27,9 +27,6 @@ tests/test_strides.py::test_strides_1arg[(10,)-None-radians] tests/test_umath.py::test_umaths[('floor_divide', 'ff')] -tests/third_party/cupy/linalg_tests/test_eigenvalue.py::TestEigenvalue_param_0_{UPLO='U'}::test_eigh_batched -tests/third_party/cupy/linalg_tests/test_eigenvalue.py::TestEigenvalue_param_1_{UPLO='L'}::test_eigh_batched - tests/third_party/cupy/random_tests/test_distributions.py::TestDistributionsBeta_param_6_{a_shape=(3, 2), b_shape=(3, 2), shape=(4, 3, 2)}::test_beta tests/third_party/cupy/random_tests/test_distributions.py::TestDistributionsBeta_param_7_{a_shape=(3, 2), b_shape=(3, 2), shape=(3, 2)}::test_beta tests/third_party/cupy/random_tests/test_distributions.py::TestDistributionsChisquare_param_0_{df_shape=(), shape=(4, 3, 2)}::test_chisquare diff --git a/tests/third_party/cupy/linalg_tests/test_eigenvalue.py b/tests/third_party/cupy/linalg_tests/test_eigenvalue.py index 93e860e5a9f1..966706e0a940 100644 --- a/tests/third_party/cupy/linalg_tests/test_eigenvalue.py +++ b/tests/third_party/cupy/linalg_tests/test_eigenvalue.py @@ -28,7 +28,6 @@ class TestEigenvalue(unittest.TestCase): rtol=1e-3, atol=1e-4, type_check=has_support_aspect64(), - contiguous_check=False, ) def test_eigh(self, xp, dtype): if xp == numpy and dtype == numpy.float16: @@ -65,7 +64,9 @@ def test_eigh(self, xp, dtype): return w @testing.for_all_dtypes(no_bool=True, no_float16=True, no_complex=True) - @testing.numpy_cupy_allclose(rtol=1e-3, atol=1e-4, contiguous_check=False) + @testing.numpy_cupy_allclose( + rtol=1e-3, atol=1e-4, type_check=has_support_aspect64() + ) def test_eigh_batched(self, xp, dtype): a = xp.array( [ @@ -89,7 +90,9 @@ def test_eigh_batched(self, xp, dtype): return w @testing.for_complex_dtypes() - @testing.numpy_cupy_allclose(rtol=1e-3, atol=1e-4, contiguous_check=False) + @testing.numpy_cupy_allclose( + rtol=1e-3, atol=1e-4, type_check=has_support_aspect64() + ) def test_eigh_complex_batched(self, xp, dtype): a = xp.array( [ @@ -195,7 +198,7 @@ def test_eigvalsh(self, xp, dtype): *testing.product( { "UPLO": ["U", "L"], - "shape": [()], + "shape": [(), (3,), (2, 3), (4, 0), (2, 2, 3), (0, 2, 3)], } ) ) @@ -203,7 +206,7 @@ class TestEigenvalueInvalid(unittest.TestCase): def test_eigh_shape_error(self): for xp in (numpy, cupy): a = xp.zeros(self.shape) - with pytest.raises((numpy.linalg.LinAlgError, ValueError)): + with pytest.raises(xp.linalg.LinAlgError): xp.linalg.eigh(a, self.UPLO) def test_eigvalsh_shape_error(self): From 6c922d9f90f59266a4addae30edda27b2f312cdb Mon Sep 17 00:00:00 2001 From: Vladislav Perevezentsev Date: Thu, 15 Feb 2024 16:44:36 +0100 Subject: [PATCH 05/26] Add test_eigenvalue_symm to check queue and usm_type --- tests/test_sycl_queue.py | 76 +++++++++++++++++++++++++++------------- tests/test_usm_type.py | 47 +++++++++++++++++++++++++ 2 files changed, 99 insertions(+), 24 deletions(-) diff --git a/tests/test_sycl_queue.py b/tests/test_sycl_queue.py index 479e96e02293..88d678340f2d 100644 --- a/tests/test_sycl_queue.py +++ b/tests/test_sycl_queue.py @@ -1087,43 +1087,71 @@ def test_eig(device): assert_sycl_queue_equal(dpnp_vec_queue, expected_queue) -@pytest.mark.usefixtures("allow_fall_back_on_numpy") +@pytest.mark.parametrize( + "func", + [ + "eigh", + "eigvalsh", + ], +) +@pytest.mark.parametrize( + "shape", + [ + (4, 4), + (0, 0), + (2, 3, 3), + (0, 2, 2), + (1, 0, 0), + ], + ids=[ + "(4, 4)", + "(0, 0)", + "(2, 3, 3)", + "(0, 2, 2)", + "(1, 0, 0)", + ], +) @pytest.mark.parametrize( "device", valid_devices, ids=[device.filter_string for device in valid_devices], ) -def test_eigh(device): - size = 4 +def test_eigenvalue_symm(func, shape, device): 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 + numpy.random.seed(81) + a = numpy.random.randn(*shape).astype(dtype) + if a.size > 0: + if a.ndim > 2: + for i in range(a.shape[0]): + a[i] = numpy.conj(a[i].T) @ a[i] + else: + a = numpy.conj(a.T) @ a - dpnp_val, dpnp_vec = dpnp.linalg.eigh(dpnp_data) - numpy_val, numpy_vec = numpy.linalg.eigh(numpy_data) + dp_a = dpnp.array(a, device=device) - assert_allclose(dpnp_val, numpy_val, rtol=1e-05, atol=1e-05) - assert_allclose(dpnp_vec, numpy_vec, rtol=1e-05, atol=1e-05) + expected_queue = dp_a.get_array().sycl_queue - 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 + if func == "eigh": + dp_val, dp_vec = dpnp.linalg.eigh(dp_a) + np_val, np_vec = numpy.linalg.eigh(a) - 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 + assert_allclose(dp_vec, np_vec, rtol=1e-05, atol=1e-05) + assert dp_vec.shape == np_vec.shape + + dpnp_vec_queue = dp_vec.get_array().sycl_queue + # 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) + assert_allclose(dp_val, np_val, rtol=1e-05, atol=1e-05) + assert dp_val.shape == np_val.shape + + dpnp_val_queue = dp_val.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( diff --git a/tests/test_usm_type.py b/tests/test_usm_type.py index 21dfb3cde678..2ffd48084dbc 100644 --- a/tests/test_usm_type.py +++ b/tests/test_usm_type.py @@ -650,6 +650,53 @@ def test_clip(usm_type): assert x.usm_type == y.usm_type +@pytest.mark.parametrize( + "func", + [ + "eigh", + "eigvalsh", + ], +) +@pytest.mark.parametrize( + "shape", + [ + (4, 4), + (0, 0), + (2, 3, 3), + (0, 2, 2), + (1, 0, 0), + ], + ids=[ + "(4, 4)", + "(0, 0)", + "(2, 3, 3)", + "(0, 2, 2)", + "(1, 0, 0)", + ], +) +@pytest.mark.parametrize("usm_type", list_of_usm_types, ids=list_of_usm_types) +def test_eigenvalue_symm(func, shape, usm_type): + dp.random.seed(81) + a = dp.random.randn(*shape) + if a.size > 0: + if a.ndim > 2: + for i in range(a.shape[0]): + a[i] = dp.conj(a[i].T) @ a[i] + else: + a = dp.conj(a.T) @ a + + a = dp.array(a, usm_type=usm_type) + + if func == "eigh": + dp_val, dp_vec = dp.linalg.eigh(a) + assert a.usm_type == dp_vec.usm_type + + else: # eighvalsh + dp_val = dp.linalg.eigvalsh(a) + + assert a.usm_type == dp_val.usm_type + + @pytest.mark.parametrize( "usm_type_matrix", list_of_usm_types, ids=list_of_usm_types ) From 2f01aed42531b10da77569b8eb1428bc0fd54e4c Mon Sep 17 00:00:00 2001 From: Vladislav Perevezentsev Date: Thu, 15 Feb 2024 17:10:03 +0100 Subject: [PATCH 06/26] Update test_linalg.py --- tests/test_linalg.py | 57 ++++++++++++++++---------------------------- 1 file changed, 21 insertions(+), 36 deletions(-) diff --git a/tests/test_linalg.py b/tests/test_linalg.py index 8e32b867b851..92ef754f2f93 100644 --- a/tests/test_linalg.py +++ b/tests/test_linalg.py @@ -383,47 +383,32 @@ def test_eig_arange(type, size): assert_allclose(dpnp_vec, np_vec, rtol=1e-05, atol=1e-05) -@pytest.mark.parametrize("type", get_all_dtypes(no_bool=True, no_none=True)) -@pytest.mark.parametrize("size", [2, 4, 8]) -def test_eigh_arange(type, size): - if dpctl.get_current_device_type() != dpctl.device_type.gpu: - pytest.skip( - "eig function doesn't work on CPU: https://github.com/IntelPython/dpnp/issues/1005" - ) - 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)) +class TestEigenvalueSymm: + @pytest.mark.parametrize( + "func", + [ + "eigh", + "eigvalsh", + ], ) - symm = symm_orig - dpnp_symm_orig = inp.array(symm) - dpnp_symm = dpnp_symm_orig - - dpnp_val, dpnp_vec = inp.linalg.eigh(dpnp_symm) - np_val, np_vec = numpy.linalg.eigh(symm) - - # DPNP sort val/vec by abs value - vvsort(dpnp_val, dpnp_vec, size, inp) + def test_eigenvalue_errors(self, func): + a_dp = inp.array([[1, 3], [3, 2]], dtype="float32") - # NP sort val/vec by abs value - vvsort(np_val, np_vec, size, numpy) - - # NP change sign of vectors - for i in range(np_vec.shape[1]): - if (np_vec[0, i] * dpnp_vec[0, i]).asnumpy() < 0: - np_vec[:, i] = -np_vec[:, i] + # unsupported type + a_np = inp.asnumpy(a_dp) + dpnp_func = getattr(inp.linalg, func) + assert_raises(TypeError, dpnp_func, a_np) - assert_array_equal(symm_orig, symm) - assert_array_equal(dpnp_symm_orig, dpnp_symm) + # a.ndim < 2 + a_dp_ndim_1 = a_dp.flatten() + assert_raises(inp.linalg.LinAlgError, dpnp_func, a_dp_ndim_1) - 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 + # a is not square + a_dp = inp.ones((2, 3)) + assert_raises(inp.linalg.LinAlgError, dpnp_func, a_dp) - assert_allclose(dpnp_val, np_val, rtol=1e-05, atol=1e-04) - assert_allclose(dpnp_vec, np_vec, rtol=1e-05, atol=1e-04) + # invalid UPLO + assert_raises(ValueError, dpnp_func, a_dp, "N") @pytest.mark.parametrize("type", get_all_dtypes(no_bool=True, no_complex=True)) From 3279e60b3cfd712d35ee2fcc7d2a384e3f385673 Mon Sep 17 00:00:00 2001 From: Vladislav Perevezentsev Date: Sun, 18 Feb 2024 14:20:52 +0100 Subject: [PATCH 07/26] Optimize dpnp_eigh logic for eigen_mode='N' --- dpnp/linalg/dpnp_utils_linalg.py | 77 ++++++++++++++++++++++++-------- 1 file changed, 59 insertions(+), 18 deletions(-) diff --git a/dpnp/linalg/dpnp_utils_linalg.py b/dpnp/linalg/dpnp_utils_linalg.py index 676d0bc886ab..e80de6b8a5c0 100644 --- a/dpnp/linalg/dpnp_utils_linalg.py +++ b/dpnp/linalg/dpnp_utils_linalg.py @@ -736,8 +736,10 @@ def dpnp_eigh(a, UPLO, eigen_mode="V"): """ dpnp_eigh(a, UPLO, eigen_mode="V") - Return the eigenvalues and eigenvectors of a complex Hermitian + Compute the eigenvalues and eigenvectors of a complex Hermitian (conjugate symmetric) or a real symmetric matrix. + Can compute both eigenvalues and eigenvectors (`eigen_mode="V"`) or + only eigenvalues (`eigen_mode="N"`). The main calculation is done by calling an extension function for LAPACK library of OneMKL. Depending on input type of `a` array, @@ -776,9 +778,34 @@ def dpnp_eigh(a, UPLO, eigen_mode="V"): # need to loop over the 1st dimension to get eigenvalues and eigenvectors of 3d matrix A batch_size = a.shape[0] if batch_size == 0: - return w, dpnp.empty_like(a, dtype=v_type) + return ( + (w, dpnp.empty_like(a, dtype=v_type)) + if eigen_mode == "V" + else w + ) eig_vecs = [None] * batch_size + + # When `eigen_mode == "N"` (jobz == 0), OneMKL LAPACK does not overwrite the input array. + # If the input array 'a' is already F-contiguous and matches the target data type, + # we avoid unnecessary memory allocation and data copying. + if eigen_mode == "N" and a_order == "F" and a.dtype == v_type: + ht_list_ev = [None] * batch_size + for i in range(batch_size): + # call LAPACK extension function to get eigenvalues of a portion of matrix A + ht_list_ev[i], _ = getattr(li, lapack_func)( + a_sycl_queue, + jobz, + uplo, + a[i].get_array(), + w[i].get_array(), + depends=[], + ) + + ht_list_ev.wait() + + return w + ht_list_ev = [None] * batch_size * 2 for i in range(batch_size): # oneMKL LAPACK assumes fortran-like array as input, so @@ -804,18 +831,32 @@ def dpnp_eigh(a, UPLO, eigen_mode="V"): dpctl.SyclEvent.wait_for(ht_list_ev) - # combine the list of eigenvectors into a single array - v = dpnp.array(eig_vecs, order=a_order) - return w, v + if eigen_mode == "V": + # combine the list of eigenvectors into a single array + v = dpnp.array(eig_vecs, order=a_order) + return w, v + return w + else: - # oneMKL LAPACK assumes fortran-like array as input, so - # allocate a memory with 'F' order for dpnp array of eigenvectors - v = dpnp.empty_like(a, order="F", dtype=v_type) + ht_list_ev = [] + copy_ev = None - # use DPCTL tensor function to fill the array of eigenvectors with content of input array - ht_copy_ev, copy_ev = ti._copy_usm_ndarray_into_usm_ndarray( - src=a_usm_arr, dst=v.get_array(), sycl_queue=a_sycl_queue - ) + # When `eigen_mode == "N"` (jobz == 0), OneMKL LAPACK does not overwrite the input array. + # If the input array 'a' is already F-contiguous and matches the target data type, + # we avoid unnecessary memory allocation and data copying. + if eigen_mode == "N" and a_order == "F" and a.dtype == v_type: + v = a + + else: + # oneMKL LAPACK assumes fortran-like array as input, so + # allocate a memory with 'F' order for dpnp array of eigenvectors + v = dpnp.empty_like(a, order="F", dtype=v_type) + + # use DPCTL tensor function to fill the array of eigenvectors with content of input array + ht_copy_ev, copy_ev = ti._copy_usm_ndarray_into_usm_ndarray( + src=a_usm_arr, dst=v.get_array(), sycl_queue=a_sycl_queue + ) + ht_list_ev.append(ht_copy_ev) # allocate a memory for dpnp array of eigenvalues w = dpnp.empty_like( @@ -833,8 +874,9 @@ def dpnp_eigh(a, UPLO, eigen_mode="V"): w.get_array(), depends=[copy_ev], ) + ht_list_ev.append(ht_lapack_ev) - if a_order != "F": + if eigen_mode == "V" and a_order != "F": # need to align order of eigenvectors with one of input matrix A out_v = dpnp.empty_like(v, order=a_order) ht_copy_out_ev, _ = ti._copy_usm_ndarray_into_usm_ndarray( @@ -843,14 +885,13 @@ def dpnp_eigh(a, UPLO, eigen_mode="V"): sycl_queue=a_sycl_queue, depends=[lapack_ev], ) - ht_copy_out_ev.wait() + ht_list_ev.append(ht_copy_out_ev) else: out_v = v - ht_lapack_ev.wait() - ht_copy_ev.wait() + dpctl.SyclEvent.wait_for(ht_list_ev) - return w, out_v + return (w, out_v) if eigen_mode == "V" else w def dpnp_eigvalsh(a, UPLO): @@ -865,7 +906,7 @@ def dpnp_eigvalsh(a, UPLO): res_type = _real_type(_common_type(a)) return dpnp.empty_like(a, shape=a.shape[:-1], dtype=res_type) - return dpnp_eigh(a, UPLO=UPLO, eigen_mode="N")[0] + return dpnp_eigh(a, UPLO=UPLO, eigen_mode="N") def dpnp_inv_batched(a, res_type): From 40e3f1999e28a388fb902b707b393ab80eb5ab8e Mon Sep 17 00:00:00 2001 From: Vladislav Perevezentsev Date: Sun, 18 Feb 2024 14:42:03 +0100 Subject: [PATCH 08/26] Add a logic with a.size==0 to dpnp_eigh --- dpnp/linalg/dpnp_utils_linalg.py | 23 +++++++++++++---------- 1 file changed, 13 insertions(+), 10 deletions(-) diff --git a/dpnp/linalg/dpnp_utils_linalg.py b/dpnp/linalg/dpnp_utils_linalg.py index e80de6b8a5c0..c6829fe4e1a6 100644 --- a/dpnp/linalg/dpnp_utils_linalg.py +++ b/dpnp/linalg/dpnp_utils_linalg.py @@ -736,9 +736,9 @@ def dpnp_eigh(a, UPLO, eigen_mode="V"): """ dpnp_eigh(a, UPLO, eigen_mode="V") - Compute the eigenvalues and eigenvectors of a complex Hermitian + Return the eigenvalues and eigenvectors of a complex Hermitian (conjugate symmetric) or a real symmetric matrix. - Can compute both eigenvalues and eigenvectors (`eigen_mode="V"`) or + Can return both eigenvalues and eigenvectors (`eigen_mode="V"`) or only eigenvalues (`eigen_mode="N"`). The main calculation is done by calling an extension function @@ -747,6 +747,17 @@ def dpnp_eigh(a, UPLO, eigen_mode="V"): """ + # get resulting type of arrays with eigenvalues and eigenvectors + v_type = _common_type(a) + w_type = _real_type(v_type) + + if a.size == 0: + w = dpnp.empty_like(a, shape=a.shape[:-1], dtype=w_type) + if eigen_mode == "V": + v = dpnp.empty_like(a, dtype=v_type) + return w, v + return w + a_sycl_queue = a.sycl_queue a_order = "C" if a.flags.c_contiguous else "F" a_usm_arr = dpnp.get_usm_ndarray(a) @@ -758,10 +769,6 @@ def dpnp_eigh(a, UPLO, eigen_mode="V"): jobz = _jobz[eigen_mode] uplo = _upper_lower[UPLO] - # get resulting type of arrays with eigenvalues and eigenvectors - v_type = _common_type(a) - w_type = _real_type(v_type) - # Get LAPACK function (_syevd for real or _heevd for complex data types) # to compute all eigenvalues and, optionally, all eigenvectors lapack_func = ( @@ -902,10 +909,6 @@ def dpnp_eigvalsh(a, UPLO): """ - if a.size == 0: - res_type = _real_type(_common_type(a)) - return dpnp.empty_like(a, shape=a.shape[:-1], dtype=res_type) - return dpnp_eigh(a, UPLO=UPLO, eigen_mode="N") From ed621853029e6c93bd84b0922d0da96faa2eefa2 Mon Sep 17 00:00:00 2001 From: Vladislav Perevezentsev Date: Sun, 18 Feb 2024 14:48:04 +0100 Subject: [PATCH 09/26] Update docstrings for eigh and eigvalsh --- dpnp/linalg/dpnp_iface_linalg.py | 14 ++++++++------ 1 file changed, 8 insertions(+), 6 deletions(-) diff --git a/dpnp/linalg/dpnp_iface_linalg.py b/dpnp/linalg/dpnp_iface_linalg.py index 945d03df4c6f..af2e97fb91d6 100644 --- a/dpnp/linalg/dpnp_iface_linalg.py +++ b/dpnp/linalg/dpnp_iface_linalg.py @@ -265,9 +265,11 @@ def eigh(a, UPLO="L"): See Also -------- - :obj:`dpnp.linalg.eig` : Return the eigenvalues and right eigenvectors of - a general matrix. - :obj:`dpnp.linalg.eigvals` : Return the eigenvalues of a general matrix. + :obj:`dpnp.linalg.eigvalsh` : Compute the eigenvalues of a complex Hermitian or + real symmetric matrix. + :obj:`dpnp.linalg.eig` : Compute the eigenvalues and right eigenvectors of + a square array. + :obj:`dpnp.linalg.eigvals` : Compute the eigenvalues of a general matrix. Examples -------- @@ -326,7 +328,7 @@ def eigvalsh(a, UPLO="L"): """ eigvalsh(a, UPLO="L") - Return the eigenvalues of a complex Hermitian or real symmetric matrix. + Compute the eigenvalues of a complex Hermitian or real symmetric matrix. Main difference from :obj:`dpnp.linalg.eigh`: the eigenvectors are not computed. @@ -355,8 +357,8 @@ def eigvalsh(a, UPLO="L"): -------- :obj:`dpnp.linalg.eigh` : Return the eigenvalues and eigenvectors of a complex Hermitian (conjugate symmetric) or a real symmetric matrix. - :obj:`dpnp.linalg.eigvals` : Return the eigenvalues of a general matrix. - :obj:`dpnp.linalg.eig` : Return the eigenvalues and right eigenvectors of + :obj:`dpnp.linalg.eigvals` : Compute the eigenvalues of a general matrix. + :obj:`dpnp.linalg.eig` : Compute the eigenvalues and right eigenvectors of a general matrix. Examples From ed8f688ef4052f022c4289c61aab76acab3f4761 Mon Sep 17 00:00:00 2001 From: Vladislav Perevezentsev Date: Sun, 18 Feb 2024 14:54:12 +0100 Subject: [PATCH 10/26] Update cupy tests --- .../cupy/linalg_tests/test_eigenvalue.py | 20 ++++++------------- 1 file changed, 6 insertions(+), 14 deletions(-) diff --git a/tests/third_party/cupy/linalg_tests/test_eigenvalue.py b/tests/third_party/cupy/linalg_tests/test_eigenvalue.py index 966706e0a940..18107d945d94 100644 --- a/tests/third_party/cupy/linalg_tests/test_eigenvalue.py +++ b/tests/third_party/cupy/linalg_tests/test_eigenvalue.py @@ -1,5 +1,3 @@ -import unittest - import numpy import pytest @@ -22,7 +20,7 @@ def _get_hermitian(xp, a, UPLO): } ) ) -class TestEigenvalue(unittest.TestCase): +class TestEigenvalue: @testing.for_all_dtypes() @testing.numpy_cupy_allclose( rtol=1e-3, @@ -90,9 +88,7 @@ def test_eigh_batched(self, xp, dtype): return w @testing.for_complex_dtypes() - @testing.numpy_cupy_allclose( - rtol=1e-3, atol=1e-4, type_check=has_support_aspect64() - ) + @testing.numpy_cupy_allclose(rtol=1e-3, atol=1e-4) def test_eigh_complex_batched(self, xp, dtype): a = xp.array( [ @@ -145,9 +141,7 @@ def test_eigvalsh_batched(self, xp, dtype): return w @testing.for_complex_dtypes() - @testing.numpy_cupy_allclose( - rtol=1e-3, atol=1e-4, type_check=has_support_aspect64() - ) + @testing.numpy_cupy_allclose(rtol=1e-3, atol=1e-4) def test_eigvalsh_complex(self, xp, dtype): a = xp.array([[1, 2j, 3], [4j, 5, 6j], [7, 8j, 9]], dtype) w = xp.linalg.eigvalsh(a, UPLO=self.UPLO) @@ -156,9 +150,7 @@ def test_eigvalsh_complex(self, xp, dtype): return w @testing.for_complex_dtypes() - @testing.numpy_cupy_allclose( - rtol=1e-3, atol=1e-4, type_check=has_support_aspect64() - ) + @testing.numpy_cupy_allclose(rtol=1e-3, atol=1e-4) def test_eigvalsh_complex_batched(self, xp, dtype): a = xp.array( [ @@ -178,7 +170,7 @@ def test_eigvalsh_complex_batched(self, xp, dtype): {"UPLO": ["U", "L"], "shape": [(0, 0), (2, 0, 0), (0, 3, 3)]} ) ) -class TestEigenvalueEmpty(unittest.TestCase): +class TestEigenvalueEmpty: @testing.for_dtypes("ifdFD") @testing.numpy_cupy_allclose(type_check=has_support_aspect64()) def test_eigh(self, xp, dtype): @@ -202,7 +194,7 @@ def test_eigvalsh(self, xp, dtype): } ) ) -class TestEigenvalueInvalid(unittest.TestCase): +class TestEigenvalueInvalid: def test_eigh_shape_error(self): for xp in (numpy, cupy): a = xp.zeros(self.shape) From 562dce4c43c859c5c6740a8b8fb67788c1ca8e7e Mon Sep 17 00:00:00 2001 From: Vladislav Perevezentsev Date: Sun, 18 Feb 2024 15:02:13 +0100 Subject: [PATCH 11/26] Add support for both registers to UPLO --- dpnp/linalg/dpnp_iface_linalg.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/dpnp/linalg/dpnp_iface_linalg.py b/dpnp/linalg/dpnp_iface_linalg.py index af2e97fb91d6..4ccffba55298 100644 --- a/dpnp/linalg/dpnp_iface_linalg.py +++ b/dpnp/linalg/dpnp_iface_linalg.py @@ -290,6 +290,7 @@ def eigh(a, UPLO="L"): check_stacked_2d(a) check_stacked_square(a) + UPLO = UPLO.upper() if UPLO not in ("L", "U"): raise ValueError("UPLO argument must be 'L' or 'U'") @@ -375,6 +376,7 @@ def eigvalsh(a, UPLO="L"): check_stacked_2d(a) check_stacked_square(a) + UPLO = UPLO.upper() if UPLO not in ("L", "U"): raise ValueError("UPLO argument must be 'L' or 'U'") From 015eb4ae3b4040cd76717c58492499a261269949 Mon Sep 17 00:00:00 2001 From: Vladislav Perevezentsev Date: Sun, 18 Feb 2024 15:05:04 +0100 Subject: [PATCH 12/26] Call dpnp_eigh instead of dpnp_eigvalsh --- dpnp/linalg/dpnp_iface_linalg.py | 3 +-- dpnp/linalg/dpnp_utils_linalg.py | 12 ------------ 2 files changed, 1 insertion(+), 14 deletions(-) diff --git a/dpnp/linalg/dpnp_iface_linalg.py b/dpnp/linalg/dpnp_iface_linalg.py index 4ccffba55298..cea82ea93724 100644 --- a/dpnp/linalg/dpnp_iface_linalg.py +++ b/dpnp/linalg/dpnp_iface_linalg.py @@ -50,7 +50,6 @@ dpnp_cholesky, dpnp_det, dpnp_eigh, - dpnp_eigvalsh, dpnp_inv, dpnp_qr, dpnp_slogdet, @@ -380,7 +379,7 @@ def eigvalsh(a, UPLO="L"): if UPLO not in ("L", "U"): raise ValueError("UPLO argument must be 'L' or 'U'") - return dpnp_eigvalsh(a, UPLO=UPLO) + return dpnp_eigh(a, UPLO=UPLO, eigen_mode="N") def inv(a): diff --git a/dpnp/linalg/dpnp_utils_linalg.py b/dpnp/linalg/dpnp_utils_linalg.py index c6829fe4e1a6..1c14fa88ced9 100644 --- a/dpnp/linalg/dpnp_utils_linalg.py +++ b/dpnp/linalg/dpnp_utils_linalg.py @@ -38,7 +38,6 @@ "dpnp_cholesky", "dpnp_det", "dpnp_eigh", - "dpnp_eigvalsh", "dpnp_inv", "dpnp_qr", "dpnp_slogdet", @@ -901,17 +900,6 @@ def dpnp_eigh(a, UPLO, eigen_mode="V"): return (w, out_v) if eigen_mode == "V" else w -def dpnp_eigvalsh(a, UPLO): - """ - dpnp_eigvalsh(a, UPLO) - - Return the eigenvalues of a complex Hermitian or real symmetric matrix. - - """ - - return dpnp_eigh(a, UPLO=UPLO, eigen_mode="N") - - def dpnp_inv_batched(a, res_type): """ dpnp_inv_batched(a, res_type) From 17691c4cb873e3fdb3bf6447b45d0cb9965b9d44 Mon Sep 17 00:00:00 2001 From: Vladislav Perevezentsev Date: Sun, 18 Feb 2024 18:08:07 +0100 Subject: [PATCH 13/26] Add get_symm_herm_numpy_array func to helper.py --- tests/helper.py | 25 +++++++++++++++++++++++++ tests/test_sycl_queue.py | 21 ++++++++++----------- tests/test_usm_type.py | 16 ++++------------ 3 files changed, 39 insertions(+), 23 deletions(-) diff --git a/tests/helper.py b/tests/helper.py index 2a2873afdcea..b10ba7524f42 100644 --- a/tests/helper.py +++ b/tests/helper.py @@ -166,6 +166,31 @@ def get_all_dtypes( return dtypes +def get_symm_herm_numpy_array(shape, dtype=None): + """ + Generates a real symmetric or a complex Hermitian numpy array of + the specified shape and data type. + + Note: + For arrays with more than 2 dimensions, it ensures symmetry(or Hermitian property + for complex data type) for each sub-array. + + """ + + numpy.random.seed(81) + a = numpy.random.randn(*shape).astype(dtype) + if numpy.issubdtype(dtype, numpy.complexfloating): + a += 1j * numpy.random.randn(*shape) + + if a.size > 0: + if a.ndim > 2: + for i in range(a.shape[0]): + a[i] = numpy.conj(a[i].T) @ a[i] + else: + a = numpy.conj(a.T) @ a + return a + + def is_cpu_device(device=None): """ Return True if a test is running on CPU device, False otherwise. diff --git a/tests/test_sycl_queue.py b/tests/test_sycl_queue.py index 88d678340f2d..3f44b6cc7e1a 100644 --- a/tests/test_sycl_queue.py +++ b/tests/test_sycl_queue.py @@ -8,7 +8,12 @@ import dpnp from dpnp.dpnp_array import dpnp_array -from .helper import assert_dtype_allclose, get_all_dtypes, is_win_platform +from .helper import ( + assert_dtype_allclose, + get_all_dtypes, + get_symm_herm_numpy_array, + is_win_platform, +) list_of_backend_str = [ "host", @@ -1116,17 +1121,9 @@ def test_eig(device): valid_devices, ids=[device.filter_string for device in valid_devices], ) -def test_eigenvalue_symm(func, shape, device): +def test_eigenvalue(func, shape, device): dtype = dpnp.default_float_type(device) - numpy.random.seed(81) - a = numpy.random.randn(*shape).astype(dtype) - if a.size > 0: - if a.ndim > 2: - for i in range(a.shape[0]): - a[i] = numpy.conj(a[i].T) @ a[i] - else: - a = numpy.conj(a.T) @ a - + a = get_symm_herm_numpy_array(shape, dtype) dp_a = dpnp.array(a, device=device) expected_queue = dp_a.get_array().sycl_queue @@ -1137,6 +1134,7 @@ def test_eigenvalue_symm(func, shape, device): assert_allclose(dp_vec, np_vec, rtol=1e-05, atol=1e-05) assert dp_vec.shape == np_vec.shape + assert dp_vec.dtype == np_vec.dtype dpnp_vec_queue = dp_vec.get_array().sycl_queue # compare queue and device @@ -1148,6 +1146,7 @@ def test_eigenvalue_symm(func, shape, device): assert_allclose(dp_val, np_val, rtol=1e-05, atol=1e-05) assert dp_val.shape == np_val.shape + assert dp_val.dtype == np_val.dtype dpnp_val_queue = dp_val.get_array().sycl_queue # compare queue and device diff --git a/tests/test_usm_type.py b/tests/test_usm_type.py index 2ffd48084dbc..b2509da52eea 100644 --- a/tests/test_usm_type.py +++ b/tests/test_usm_type.py @@ -8,7 +8,7 @@ import dpnp as dp -from .helper import assert_dtype_allclose +from .helper import assert_dtype_allclose, get_symm_herm_numpy_array list_of_usm_types = ["device", "shared", "host"] @@ -675,17 +675,9 @@ def test_clip(usm_type): ], ) @pytest.mark.parametrize("usm_type", list_of_usm_types, ids=list_of_usm_types) -def test_eigenvalue_symm(func, shape, usm_type): - dp.random.seed(81) - a = dp.random.randn(*shape) - if a.size > 0: - if a.ndim > 2: - for i in range(a.shape[0]): - a[i] = dp.conj(a[i].T) @ a[i] - else: - a = dp.conj(a.T) @ a - - a = dp.array(a, usm_type=usm_type) +def test_eigenvalue(func, shape, usm_type): + a_np = get_symm_herm_numpy_array(shape) + a = dp.array(a_np, usm_type=usm_type) if func == "eigh": dp_val, dp_vec = dp.linalg.eigh(a) From f8eb398da00e39b124410c591057575801796795 Mon Sep 17 00:00:00 2001 From: Vladislav Perevezentsev Date: Sun, 18 Feb 2024 19:31:05 +0100 Subject: [PATCH 14/26] Update dpnp tests --- tests/helper.py | 2 +- tests/test_linalg.py | 13 +++---------- tests/test_sycl_queue.py | 6 +++--- 3 files changed, 7 insertions(+), 14 deletions(-) diff --git a/tests/helper.py b/tests/helper.py index b10ba7524f42..bba5edf3f170 100644 --- a/tests/helper.py +++ b/tests/helper.py @@ -179,7 +179,7 @@ def get_symm_herm_numpy_array(shape, dtype=None): numpy.random.seed(81) a = numpy.random.randn(*shape).astype(dtype) - if numpy.issubdtype(dtype, numpy.complexfloating): + if numpy.issubdtype(a.dtype, numpy.complexfloating): a += 1j * numpy.random.randn(*shape) if a.size > 0: diff --git a/tests/test_linalg.py b/tests/test_linalg.py index 2cb01486727a..ed948ea39ce6 100644 --- a/tests/test_linalg.py +++ b/tests/test_linalg.py @@ -15,6 +15,7 @@ assert_dtype_allclose, get_all_dtypes, get_float_complex_dtypes, + get_symm_herm_numpy_array, has_support_aspect64, is_cpu_device, ) @@ -1126,11 +1127,7 @@ def test_svd(self, dtype, shape): ids=["(2, 2)", "(16, 16)"], ) def test_svd_hermitian(self, dtype, compute_vt, shape): - a = numpy.random.randn(*shape).astype(dtype) - if numpy.issubdtype(dtype, numpy.complexfloating): - a += 1j * numpy.random.randn(*shape) - a = (a + a.conj().T) / 2 - + a = get_symm_herm_numpy_array(shape, dtype) dp_a = inp.array(a) if compute_vt: @@ -1235,11 +1232,7 @@ def test_pinv(self, dtype, shape): ids=["(2, 2)", "(16, 16)"], ) def test_pinv_hermitian(self, dtype, shape): - a = numpy.random.randn(*shape).astype(dtype) - if numpy.issubdtype(dtype, numpy.complexfloating): - a += 1j * numpy.random.randn(*shape) - a = (a + a.conj().T) / 2 - + a = get_symm_herm_numpy_array(shape, dtype) a_dp = inp.array(a) B = numpy.linalg.pinv(a, hermitian=True) diff --git a/tests/test_sycl_queue.py b/tests/test_sycl_queue.py index 5039d7a7fbdc..857493d6f687 100644 --- a/tests/test_sycl_queue.py +++ b/tests/test_sycl_queue.py @@ -1726,11 +1726,11 @@ def test_slogdet(shape, is_empty, device): ) def test_pinv(shape, hermitian, rcond_as_array, device): numpy.random.seed(81) + dtype = dpnp.default_float_type(device) if hermitian: - a_np = numpy.random.randn(*shape) + 1j * numpy.random.randn(*shape) - a_np = numpy.conj(a_np.T) @ a_np + a_np = get_symm_herm_numpy_array(shape, dtype) else: - a_np = numpy.random.randn(*shape) + a_np = numpy.random.randn(*shape).astype(dtype) a_dp = dpnp.array(a_np, device=device) From 47aba16adb658c20982b061af3fdbc3ca85630f6 Mon Sep 17 00:00:00 2001 From: Vladislav Perevezentsev Date: Tue, 20 Feb 2024 10:40:16 +0100 Subject: [PATCH 15/26] Add generate_random_numpy_array func to helper.py --- tests/helper.py | 47 ++++++++++++++++++++++++++++++++++------ tests/test_linalg.py | 45 ++++++++++++++++---------------------- tests/test_sycl_queue.py | 12 +++++----- tests/test_usm_type.py | 13 ++++------- 4 files changed, 68 insertions(+), 49 deletions(-) diff --git a/tests/helper.py b/tests/helper.py index bba5edf3f170..c67d92196341 100644 --- a/tests/helper.py +++ b/tests/helper.py @@ -166,23 +166,56 @@ def get_all_dtypes( return dtypes -def get_symm_herm_numpy_array(shape, dtype=None): +def generate_random_numpy_array( + shape, dtype=None, hermitian=False, use_seed=False, seed_value=81 +): """ - Generates a real symmetric or a complex Hermitian numpy array of - the specified shape and data type. + Generate a random numpy array with the specified shape and dtype. + + If required, the array can be made Hermitian (for complex data types) or + symmetric (for real data types). + + Parameters + ---------- + shape : tuple + Shape of the generated array. + dtype : str or dtype, optional + Desired data-type for the output array. + If not specified, data type will be determined by numpy. + Default : None + hermitian : bool, optional + If True, generates a Hermitian (symmetric if `dtype` is real) matrix. + Default : False + set_seed : bool, optional + If True, the random number generator seed will be set to `seed_value`. + Default : False. + seed_value : int, optional + The seed value to initialize the random number generator. + Only used if `set_seed` is True. + Default : `81` to avoid generating a singular matrix. + + Returns + ------- + out : numpy.ndarray + A random numpy array of the specified shape and dtype. + The array is Hermitian or symmetric if `hermitian` is True. Note: - For arrays with more than 2 dimensions, it ensures symmetry(or Hermitian property - for complex data type) for each sub-array. + For arrays with more than 2 dimensions, the Hermitian or + symmetric property is ensured for each 2D sub-array. """ - numpy.random.seed(81) + if use_seed: + numpy.random.seed(seed_value) + a = numpy.random.randn(*shape).astype(dtype) if numpy.issubdtype(a.dtype, numpy.complexfloating): + if use_seed: + numpy.random.seed(seed_value) a += 1j * numpy.random.randn(*shape) - if a.size > 0: + if hermitian and a.size > 0: if a.ndim > 2: for i in range(a.shape[0]): a[i] = numpy.conj(a[i].T) @ a[i] diff --git a/tests/test_linalg.py b/tests/test_linalg.py index ed948ea39ce6..4ed590f8a9a4 100644 --- a/tests/test_linalg.py +++ b/tests/test_linalg.py @@ -13,9 +13,9 @@ from .helper import ( assert_dtype_allclose, + generate_random_numpy_array, get_all_dtypes, get_float_complex_dtypes, - get_symm_herm_numpy_array, has_support_aspect64, is_cpu_device, ) @@ -664,11 +664,6 @@ def test_norm3(array, ord, axis): class TestQr: - # Set numpy.random.seed for test methods to prevent - # random generation of the input singular matrix - def setup_method(self): - numpy.random.seed(81) - # TODO: New packages that fix issue CMPLRLLVM-53771 are only available in internal CI. # Skip the tests on cpu until these packages are available for the external CI. # Specifically dpcpp_linux-64>=2024.1.0 @@ -693,9 +688,9 @@ def setup_method(self): ids=["r", "raw", "complete", "reduced"], ) def test_qr(self, dtype, shape, mode): - a = numpy.random.randn(*shape).astype(dtype) - if numpy.issubdtype(dtype, numpy.complexfloating): - a += 1j * numpy.random.randn(*shape) + # Set use_seed=True to prevent + # random generation of the input singular matrix + a = generate_random_numpy_array(shape, dtype, use_seed=True) ia = inp.array(a) if mode == "r": @@ -765,7 +760,7 @@ def test_qr_empty(self, dtype, shape, mode): ids=["r", "raw", "complete", "reduced"], ) def test_qr_strides(self, mode): - a = numpy.random.randn(5, 5) + a = generate_random_numpy_array((5, 5)) ia = inp.array(a) # positive strides @@ -1025,11 +1020,6 @@ def test_slogdet_errors(self): class TestSvd: - # Set numpy.random.seed for test methods to prevent - # random generation of the input singular matrix - def setup_method(self): - numpy.random.seed(81) - def get_tol(self, dtype): tol = 1e-06 if dtype in (inp.float32, inp.complex64): @@ -1127,7 +1117,11 @@ def test_svd(self, dtype, shape): ids=["(2, 2)", "(16, 16)"], ) def test_svd_hermitian(self, dtype, compute_vt, shape): - a = get_symm_herm_numpy_array(shape, dtype) + # Set use_seed=True to prevent + # random generation of the input singular matrix + a = generate_random_numpy_array( + shape, dtype, hermitian=True, use_seed=True + ) dp_a = inp.array(a) if compute_vt: @@ -1165,11 +1159,6 @@ def test_svd_errors(self): class TestPinv: - # Set numpy.random.seed for test methods to prevent - # random generation of the input singular matrix - def setup_method(self): - numpy.random.seed(81) - def get_tol(self, dtype): tol = 1e-06 if dtype in (inp.float32, inp.complex64): @@ -1205,9 +1194,9 @@ def check_types_shapes(self, dp_B, np_B): ], ) def test_pinv(self, dtype, shape): - a = numpy.random.randn(*shape).astype(dtype) - if numpy.issubdtype(dtype, numpy.complexfloating): - a += 1j * numpy.random.randn(*shape) + # Set use_seed=True to prevent + # random generation of the input singular matrix + a = generate_random_numpy_array(shape, dtype, use_seed=True) a_dp = inp.array(a) B = numpy.linalg.pinv(a) @@ -1232,7 +1221,11 @@ def test_pinv(self, dtype, shape): ids=["(2, 2)", "(16, 16)"], ) def test_pinv_hermitian(self, dtype, shape): - a = get_symm_herm_numpy_array(shape, dtype) + # Set use_seed=True to prevent + # random generation of the input singular matrix + a = generate_random_numpy_array( + shape, dtype, hermitian=True, use_seed=True + ) a_dp = inp.array(a) B = numpy.linalg.pinv(a, hermitian=True) @@ -1268,7 +1261,7 @@ def test_pinv_empty(self, dtype, shape): assert_dtype_allclose(B_dp, B) def test_pinv_strides(self): - a = numpy.random.randn(5, 5) + a = generate_random_numpy_array((5, 5)) a_dp = inp.array(a) self.get_tol(a_dp.dtype) diff --git a/tests/test_sycl_queue.py b/tests/test_sycl_queue.py index 857493d6f687..be32194fa0e9 100644 --- a/tests/test_sycl_queue.py +++ b/tests/test_sycl_queue.py @@ -10,8 +10,8 @@ from .helper import ( assert_dtype_allclose, + generate_random_numpy_array, get_all_dtypes, - get_symm_herm_numpy_array, is_win_platform, ) @@ -1123,7 +1123,7 @@ def test_eig(device): ) def test_eigenvalue(func, shape, device): dtype = dpnp.default_float_type(device) - a = get_symm_herm_numpy_array(shape, dtype) + a = generate_random_numpy_array(shape, dtype, hermitian=True, use_seed=True) dp_a = dpnp.array(a, device=device) expected_queue = dp_a.get_array().sycl_queue @@ -1725,12 +1725,10 @@ def test_slogdet(shape, is_empty, device): ids=[device.filter_string for device in valid_devices], ) def test_pinv(shape, hermitian, rcond_as_array, device): - numpy.random.seed(81) dtype = dpnp.default_float_type(device) - if hermitian: - a_np = get_symm_herm_numpy_array(shape, dtype) - else: - a_np = numpy.random.randn(*shape).astype(dtype) + a_np = generate_random_numpy_array( + shape, dtype, hermitian=hermitian, use_seed=True + ) a_dp = dpnp.array(a_np, device=device) diff --git a/tests/test_usm_type.py b/tests/test_usm_type.py index f009c9b59798..42ced58e4e61 100644 --- a/tests/test_usm_type.py +++ b/tests/test_usm_type.py @@ -8,7 +8,7 @@ import dpnp as dp -from .helper import assert_dtype_allclose, get_symm_herm_numpy_array +from .helper import assert_dtype_allclose, generate_random_numpy_array list_of_usm_types = ["device", "shared", "host"] @@ -676,7 +676,7 @@ def test_clip(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 = get_symm_herm_numpy_array(shape) + a_np = generate_random_numpy_array(shape, hermitian=True) a = dp.array(a_np, usm_type=usm_type) if func == "eigh": @@ -889,14 +889,9 @@ def test_svd(usm_type, shape, full_matrices_param, compute_uv_param): ], ) def test_pinv(shape, hermitian, usm_type): - numpy.random.seed(81) - if hermitian: - a = dp.random.randn(*shape) + 1j * dp.random.randn(*shape) - a = dp.conj(a.T) @ a - else: - a = dp.random.randn(*shape) + a_np = generate_random_numpy_array(shape, hermitian=hermitian) + a = dp.array(a_np, usm_type=usm_type) - a = dp.array(a, usm_type=usm_type) B = dp.linalg.pinv(a, hermitian=hermitian) assert a.usm_type == B.usm_type From 303bdc671e554ada57e0c148a0582a312e9189c1 Mon Sep 17 00:00:00 2001 From: Vladislav Perevezentsev Date: Tue, 20 Feb 2024 16:07:40 +0100 Subject: [PATCH 16/26] Remove use_seed parameter from generate_random_numpy_array --- tests/helper.py | 14 ++++---------- tests/test_linalg.py | 16 ++++++++-------- tests/test_sycl_queue.py | 9 ++++++--- 3 files changed, 18 insertions(+), 21 deletions(-) diff --git a/tests/helper.py b/tests/helper.py index c67d92196341..f9fee1564b3a 100644 --- a/tests/helper.py +++ b/tests/helper.py @@ -167,7 +167,7 @@ def get_all_dtypes( def generate_random_numpy_array( - shape, dtype=None, hermitian=False, use_seed=False, seed_value=81 + shape, dtype=None, hermitian=False, seed_value=None ): """ Generate a random numpy array with the specified shape and dtype. @@ -186,13 +186,9 @@ def generate_random_numpy_array( hermitian : bool, optional If True, generates a Hermitian (symmetric if `dtype` is real) matrix. Default : False - set_seed : bool, optional - If True, the random number generator seed will be set to `seed_value`. - Default : False. seed_value : int, optional The seed value to initialize the random number generator. - Only used if `set_seed` is True. - Default : `81` to avoid generating a singular matrix. + Default : None Returns ------- @@ -206,13 +202,11 @@ def generate_random_numpy_array( """ - if use_seed: - numpy.random.seed(seed_value) + numpy.random.seed(seed_value) a = numpy.random.randn(*shape).astype(dtype) if numpy.issubdtype(a.dtype, numpy.complexfloating): - if use_seed: - numpy.random.seed(seed_value) + numpy.random.seed(seed_value) a += 1j * numpy.random.randn(*shape) if hermitian and a.size > 0: diff --git a/tests/test_linalg.py b/tests/test_linalg.py index 7d9c2d9a4f61..814c3fd05d49 100644 --- a/tests/test_linalg.py +++ b/tests/test_linalg.py @@ -765,9 +765,9 @@ class TestQr: ids=["r", "raw", "complete", "reduced"], ) def test_qr(self, dtype, shape, mode): - # Set use_seed=True to prevent + # Set seed_value=81 to prevent # random generation of the input singular matrix - a = generate_random_numpy_array(shape, dtype, use_seed=True) + a = generate_random_numpy_array(shape, dtype, seed_value=81) ia = inp.array(a) if mode == "r": @@ -1194,10 +1194,10 @@ def test_svd(self, dtype, shape): ids=["(2, 2)", "(16, 16)"], ) def test_svd_hermitian(self, dtype, compute_vt, shape): - # Set use_seed=True to prevent + # Set seed_value=81 to prevent # random generation of the input singular matrix a = generate_random_numpy_array( - shape, dtype, hermitian=True, use_seed=True + shape, dtype, hermitian=True, seed_value=81 ) dp_a = inp.array(a) @@ -1271,9 +1271,9 @@ def check_types_shapes(self, dp_B, np_B): ], ) def test_pinv(self, dtype, shape): - # Set use_seed=True to prevent + # Set seed_value=81 to prevent # random generation of the input singular matrix - a = generate_random_numpy_array(shape, dtype, use_seed=True) + a = generate_random_numpy_array(shape, dtype, seed_value=81) a_dp = inp.array(a) B = numpy.linalg.pinv(a) @@ -1298,10 +1298,10 @@ def test_pinv(self, dtype, shape): ids=["(2, 2)", "(16, 16)"], ) def test_pinv_hermitian(self, dtype, shape): - # Set use_seed=True to prevent + # Set seed_value=81 to prevent # random generation of the input singular matrix a = generate_random_numpy_array( - shape, dtype, hermitian=True, use_seed=True + shape, dtype, hermitian=True, seed_value=81 ) a_dp = inp.array(a) diff --git a/tests/test_sycl_queue.py b/tests/test_sycl_queue.py index 87e825a35d47..75e9efd19859 100644 --- a/tests/test_sycl_queue.py +++ b/tests/test_sycl_queue.py @@ -1123,7 +1123,9 @@ def test_eig(device): ) def test_eigenvalue(func, shape, device): dtype = dpnp.default_float_type(device) - a = generate_random_numpy_array(shape, dtype, hermitian=True, use_seed=True) + # 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) dp_a = dpnp.array(a, device=device) expected_queue = dp_a.get_array().sycl_queue @@ -1742,10 +1744,11 @@ def test_slogdet(shape, is_empty, device): ) def test_pinv(shape, hermitian, rcond_as_array, device): dtype = dpnp.default_float_type(device) + # Set seed_value=81 to prevent + # random generation of the input singular matrix a_np = generate_random_numpy_array( - shape, dtype, hermitian=hermitian, use_seed=True + shape, dtype, hermitian=hermitian, seed_value=81 ) - a_dp = dpnp.array(a_np, device=device) if rcond_as_array: From 768639a4a9cb6f17f138f62e8ace3e6392ba2f2e Mon Sep 17 00:00:00 2001 From: Vladislav Perevezentsev Date: Tue, 20 Feb 2024 16:13:04 +0100 Subject: [PATCH 17/26] Address remarks --- dpnp/linalg/dpnp_utils_linalg.py | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/dpnp/linalg/dpnp_utils_linalg.py b/dpnp/linalg/dpnp_utils_linalg.py index 7cca4910ef48..f475dcdd3de6 100644 --- a/dpnp/linalg/dpnp_utils_linalg.py +++ b/dpnp/linalg/dpnp_utils_linalg.py @@ -792,11 +792,9 @@ def dpnp_eigh(a, UPLO, eigen_mode="V"): else w ) - eig_vecs = [None] * batch_size - # When `eigen_mode == "N"` (jobz == 0), OneMKL LAPACK does not overwrite the input array. # If the input array 'a' is already F-contiguous and matches the target data type, - # we avoid unnecessary memory allocation and data copying. + # we can avoid unnecessary memory allocation and data copying. if eigen_mode == "N" and a_order == "F" and a.dtype == v_type: ht_list_ev = [None] * batch_size for i in range(batch_size): @@ -810,10 +808,11 @@ def dpnp_eigh(a, UPLO, eigen_mode="V"): depends=[], ) - ht_list_ev.wait() + dpctl.SyclEvent.wait_for(ht_list_ev) return w + eig_vecs = [None] * batch_size ht_list_ev = [None] * batch_size * 2 for i in range(batch_size): # oneMKL LAPACK assumes fortran-like array as input, so @@ -851,7 +850,7 @@ def dpnp_eigh(a, UPLO, eigen_mode="V"): # When `eigen_mode == "N"` (jobz == 0), OneMKL LAPACK does not overwrite the input array. # If the input array 'a' is already F-contiguous and matches the target data type, - # we avoid unnecessary memory allocation and data copying. + # we can avoid unnecessary memory allocation and data copying. if eigen_mode == "N" and a_order == "F" and a.dtype == v_type: v = a From 7b7bea937048b1461f61392d0f3dfe63969a79ae Mon Sep 17 00:00:00 2001 From: Vladislav Perevezentsev Date: Wed, 21 Feb 2024 11:41:17 +0100 Subject: [PATCH 18/26] Generate 4d and more array --- tests/helper.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/tests/helper.py b/tests/helper.py index f9fee1564b3a..ba6e32a247c8 100644 --- a/tests/helper.py +++ b/tests/helper.py @@ -211,8 +211,12 @@ def generate_random_numpy_array( if hermitian and a.size > 0: if a.ndim > 2: + orig_shape = a.shape + # get 3D array + a = a.reshape(-1, orig_shape[-2], orig_shape[-1]) for i in range(a.shape[0]): a[i] = numpy.conj(a[i].T) @ a[i] + a = a.reshape(orig_shape) else: a = numpy.conj(a.T) @ a return a From 6eab947570d8f1f25eb662ce1d4026a6cc6bd0c4 Mon Sep 17 00:00:00 2001 From: Vladislav Perevezentsev Date: Thu, 22 Feb 2024 10:41:37 +0100 Subject: [PATCH 19/26] Support 4d and more array for dpnp_eigh --- dpnp/linalg/dpnp_utils_linalg.py | 52 +++++++++++--------------------- 1 file changed, 18 insertions(+), 34 deletions(-) diff --git a/dpnp/linalg/dpnp_utils_linalg.py b/dpnp/linalg/dpnp_utils_linalg.py index f475dcdd3de6..7a7bf24432df 100644 --- a/dpnp/linalg/dpnp_utils_linalg.py +++ b/dpnp/linalg/dpnp_utils_linalg.py @@ -759,10 +759,6 @@ def dpnp_eigh(a, UPLO, eigen_mode="V"): return w, v return w - a_sycl_queue = a.sycl_queue - a_order = "C" if a.flags.c_contiguous else "F" - a_usm_arr = dpnp.get_usm_ndarray(a) - # `eigen_mode` can be either "N" or "V", specifying the computation mode # for OneMKL LAPACK `syevd` and `heevd` routines. # "V" (default) means both eigenvectors and eigenvalues will be calculated @@ -776,42 +772,27 @@ def dpnp_eigh(a, UPLO, eigen_mode="V"): "_heevd" if dpnp.issubdtype(v_type, dpnp.complexfloating) else "_syevd" ) + a_sycl_queue = a.sycl_queue + a_order = "C" if a.flags.c_contiguous else "F" + if a.ndim > 2: + orig_shape = a.shape + # get 3d input array by reshape + a = a.reshape(-1, orig_shape[-2], orig_shape[-1]) + a_usm_arr = dpnp.get_usm_ndarray(a) + + # allocate a memory for dpnp array of eigenvalues w = dpnp.empty_like( a, - shape=a.shape[:-1], + shape=orig_shape[:-1], dtype=w_type, ) + w_orig_shape = w.shape + # get 2d dpnp array with eigenvalues by reshape + w = w.reshape(-1, w_orig_shape[-1]) # need to loop over the 1st dimension to get eigenvalues and eigenvectors of 3d matrix A batch_size = a.shape[0] - if batch_size == 0: - return ( - (w, dpnp.empty_like(a, dtype=v_type)) - if eigen_mode == "V" - else w - ) - - # When `eigen_mode == "N"` (jobz == 0), OneMKL LAPACK does not overwrite the input array. - # If the input array 'a' is already F-contiguous and matches the target data type, - # we can avoid unnecessary memory allocation and data copying. - if eigen_mode == "N" and a_order == "F" and a.dtype == v_type: - ht_list_ev = [None] * batch_size - for i in range(batch_size): - # call LAPACK extension function to get eigenvalues of a portion of matrix A - ht_list_ev[i], _ = getattr(li, lapack_func)( - a_sycl_queue, - jobz, - uplo, - a[i].get_array(), - w[i].get_array(), - depends=[], - ) - - dpctl.SyclEvent.wait_for(ht_list_ev) - - return w - eig_vecs = [None] * batch_size ht_list_ev = [None] * batch_size * 2 for i in range(batch_size): @@ -838,15 +819,18 @@ def dpnp_eigh(a, UPLO, eigen_mode="V"): dpctl.SyclEvent.wait_for(ht_list_ev) + w = w.reshape(w_orig_shape) + if eigen_mode == "V": # combine the list of eigenvectors into a single array - v = dpnp.array(eig_vecs, order=a_order) + v = dpnp.array(eig_vecs, order=a_order).reshape(orig_shape) return w, v return w else: + a_usm_arr = dpnp.get_usm_ndarray(a) ht_list_ev = [] - copy_ev = None + copy_ev = dpctl.SyclEvent() # When `eigen_mode == "N"` (jobz == 0), OneMKL LAPACK does not overwrite the input array. # If the input array 'a' is already F-contiguous and matches the target data type, From 260b17104a2cb964be61c388f43ab5505965e51a Mon Sep 17 00:00:00 2001 From: Vladislav Perevezentsev Date: Thu, 22 Feb 2024 10:48:48 +0100 Subject: [PATCH 20/26] Update TestEigenvalue --- tests/test_linalg.py | 43 +++++++++++++++++++++++++++++++++++++++++-- 1 file changed, 41 insertions(+), 2 deletions(-) diff --git a/tests/test_linalg.py b/tests/test_linalg.py index 814c3fd05d49..e9ba31950b5f 100644 --- a/tests/test_linalg.py +++ b/tests/test_linalg.py @@ -385,7 +385,46 @@ def test_eig_arange(type, size): assert_allclose(dpnp_vec, np_vec, rtol=1e-05, atol=1e-05) -class TestEigenvalueSymm: +class TestEigenvalue: + @pytest.mark.parametrize( + "func", + [ + "eigh", + "eigvalsh", + ], + ) + @pytest.mark.parametrize( + "shape", + [(2, 2), (2, 3, 3), (2, 2, 3, 3)], + ids=["(2,2)", "(2,3,3)", "(2,2,3,3)"], + ) + @pytest.mark.parametrize("dtype", get_all_dtypes(no_bool=True)) + @pytest.mark.parametrize( + "order", + [ + "C", + "F", + ], + ) + def test_eigenvalues(self, func, shape, dtype, order): + a = generate_random_numpy_array( + shape, dtype, hermitian=True, seed_value=81 + ) + a_order = numpy.array(a, order=order) + a_dp = inp.array(a, order=order) + + if func == "eigh": + w, v = numpy.linalg.eigh(a_order) + w_dp, v_dp = inp.linalg.eigh(a_dp) + + assert_dtype_allclose(v_dp, v) + + else: # eighvalsh + w = numpy.linalg.eigvalsh(a_order) + w_dp = inp.linalg.eigvalsh(a_dp) + + assert_dtype_allclose(w_dp, w) + @pytest.mark.parametrize( "func", [ @@ -410,7 +449,7 @@ def test_eigenvalue_errors(self, func): assert_raises(inp.linalg.LinAlgError, dpnp_func, a_dp) # invalid UPLO - assert_raises(ValueError, dpnp_func, a_dp, "N") + assert_raises(ValueError, dpnp_func, a_dp, UPLO="N") @pytest.mark.parametrize("type", get_all_dtypes(no_bool=True, no_complex=True)) From 6cd5c37e4837376f9c1944fa65a08cda8b5e5942 Mon Sep 17 00:00:00 2001 From: Vladislav Perevezentsev Date: Thu, 22 Feb 2024 21:34:51 +0100 Subject: [PATCH 21/26] Fix invalid UPLO test --- tests/test_linalg.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tests/test_linalg.py b/tests/test_linalg.py index e9ba31950b5f..133d38deedca 100644 --- a/tests/test_linalg.py +++ b/tests/test_linalg.py @@ -445,8 +445,8 @@ def test_eigenvalue_errors(self, func): assert_raises(inp.linalg.LinAlgError, dpnp_func, a_dp_ndim_1) # a is not square - a_dp = inp.ones((2, 3)) - assert_raises(inp.linalg.LinAlgError, dpnp_func, a_dp) + a_dp_not_scquare = inp.ones((2, 3)) + assert_raises(inp.linalg.LinAlgError, dpnp_func, a_dp_not_scquare) # invalid UPLO assert_raises(ValueError, dpnp_func, a_dp, UPLO="N") From 08626dc72e6210c287e7719ec47a32a9a7ad6d5f Mon Sep 17 00:00:00 2001 From: Vladislav Perevezentsev Date: Mon, 25 Mar 2024 12:33:40 +0100 Subject: [PATCH 22/26] Make depends parameter optional in heevd.hpp --- dpnp/backend/extensions/lapack/heevd.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/dpnp/backend/extensions/lapack/heevd.hpp b/dpnp/backend/extensions/lapack/heevd.hpp index 89ecfe466fb8..7b3bfc05d87a 100644 --- a/dpnp/backend/extensions/lapack/heevd.hpp +++ b/dpnp/backend/extensions/lapack/heevd.hpp @@ -44,7 +44,7 @@ extern std::pair const std::int8_t upper_lower, dpctl::tensor::usm_ndarray eig_vecs, dpctl::tensor::usm_ndarray eig_vals, - const std::vector &depends); + const std::vector &depends = {}); extern void init_heevd_dispatch_table(void); } // namespace lapack From 2e0c3c4ac940854ee09f7f115e17e01327cb6367 Mon Sep 17 00:00:00 2001 From: Vladislav Perevezentsev Date: Mon, 25 Mar 2024 12:36:24 +0100 Subject: [PATCH 23/26] Add a w/a for MKLD-17201 --- dpnp/linalg/dpnp_utils_linalg.py | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/dpnp/linalg/dpnp_utils_linalg.py b/dpnp/linalg/dpnp_utils_linalg.py index 73fcfa49e4fa..bdb465b205c6 100644 --- a/dpnp/linalg/dpnp_utils_linalg.py +++ b/dpnp/linalg/dpnp_utils_linalg.py @@ -858,6 +858,7 @@ def dpnp_eigh(a, UPLO, eigen_mode="V"): a_order = "C" if a.flags.c_contiguous else "F" if a.ndim > 2: + is_cpu_device = a.sycl_device.has_aspect_cpu orig_shape = a.shape # get 3d input array by reshape a = a.reshape(-1, orig_shape[-2], orig_shape[-1]) @@ -889,6 +890,14 @@ def dpnp_eigh(a, UPLO, eigen_mode="V"): sycl_queue=a_sycl_queue, ) + # TODO: Remove this w/a when MKLD-17201 is solved. + # Waiting for a host task executing an OneMKL LAPACK syevd call + # on CPU causes deadlock due to serialization of all host tasks + # in the queue. + # We need to wait for each host tasks before calling _seyvd to avoid deadlock. + if lapack_func == "_syevd" and is_cpu_device: + ht_list_ev[2 * i].wait() + # call LAPACK extension function to get eigenvalues and eigenvectors of a portion of matrix A ht_list_ev[2 * i + 1], _ = getattr(li, lapack_func)( a_sycl_queue, From 8efb5585e3c313a2b90c273d4d07ceead30d8fff Mon Sep 17 00:00:00 2001 From: Vladislav Perevezentsev Date: Tue, 26 Mar 2024 12:18:05 +0100 Subject: [PATCH 24/26] Wait for each host_task after calling _syevd --- dpnp/linalg/dpnp_utils_linalg.py | 17 +++++++++-------- 1 file changed, 9 insertions(+), 8 deletions(-) diff --git a/dpnp/linalg/dpnp_utils_linalg.py b/dpnp/linalg/dpnp_utils_linalg.py index 05feaed36647..4a9572044336 100644 --- a/dpnp/linalg/dpnp_utils_linalg.py +++ b/dpnp/linalg/dpnp_utils_linalg.py @@ -932,14 +932,6 @@ def dpnp_eigh(a, UPLO, eigen_mode="V"): sycl_queue=a_sycl_queue, ) - # TODO: Remove this w/a when MKLD-17201 is solved. - # Waiting for a host task executing an OneMKL LAPACK syevd call - # on CPU causes deadlock due to serialization of all host tasks - # in the queue. - # We need to wait for each host tasks before calling _seyvd to avoid deadlock. - if lapack_func == "_syevd" and is_cpu_device: - ht_list_ev[2 * i].wait() - # call LAPACK extension function to get eigenvalues and eigenvectors of a portion of matrix A ht_list_ev[2 * i + 1], _ = getattr(li, lapack_func)( a_sycl_queue, @@ -950,6 +942,15 @@ def dpnp_eigh(a, UPLO, eigen_mode="V"): depends=[copy_ev], ) + # TODO: Remove this w/a when MKLD-17201 is solved. + # Waiting for a host task executing an OneMKL LAPACK syevd call + # on CPU causes deadlock due to serialization of all host tasks + # in the queue. + # We need to wait for each host tasks before calling _seyvd again + # to avoid deadlock. + if lapack_func == "_syevd" and is_cpu_device: + ht_list_ev[2 * i + 1].wait() + dpctl.SyclEvent.wait_for(ht_list_ev) w = w.reshape(w_orig_shape) From d86a957727d6dee0f3656b489f582c4d0b93e8b9 Mon Sep 17 00:00:00 2001 From: Vladislav Perevezentsev Date: Tue, 26 Mar 2024 15:57:17 +0100 Subject: [PATCH 25/26] Update test_eigenvalues in test_linalg.py --- dpnp/linalg/dpnp_utils_linalg.py | 17 ++++++++--------- tests/test_linalg.py | 26 ++++++++++++++++++++++++-- 2 files changed, 32 insertions(+), 11 deletions(-) diff --git a/dpnp/linalg/dpnp_utils_linalg.py b/dpnp/linalg/dpnp_utils_linalg.py index 4a9572044336..05feaed36647 100644 --- a/dpnp/linalg/dpnp_utils_linalg.py +++ b/dpnp/linalg/dpnp_utils_linalg.py @@ -932,6 +932,14 @@ def dpnp_eigh(a, UPLO, eigen_mode="V"): sycl_queue=a_sycl_queue, ) + # TODO: Remove this w/a when MKLD-17201 is solved. + # Waiting for a host task executing an OneMKL LAPACK syevd call + # on CPU causes deadlock due to serialization of all host tasks + # in the queue. + # We need to wait for each host tasks before calling _seyvd to avoid deadlock. + if lapack_func == "_syevd" and is_cpu_device: + ht_list_ev[2 * i].wait() + # call LAPACK extension function to get eigenvalues and eigenvectors of a portion of matrix A ht_list_ev[2 * i + 1], _ = getattr(li, lapack_func)( a_sycl_queue, @@ -942,15 +950,6 @@ def dpnp_eigh(a, UPLO, eigen_mode="V"): depends=[copy_ev], ) - # TODO: Remove this w/a when MKLD-17201 is solved. - # Waiting for a host task executing an OneMKL LAPACK syevd call - # on CPU causes deadlock due to serialization of all host tasks - # in the queue. - # We need to wait for each host tasks before calling _seyvd again - # to avoid deadlock. - if lapack_func == "_syevd" and is_cpu_device: - ht_list_ev[2 * i + 1].wait() - dpctl.SyclEvent.wait_for(ht_list_ev) w = w.reshape(w_orig_shape) diff --git a/tests/test_linalg.py b/tests/test_linalg.py index bf171a30ffa3..3977d538f09d 100644 --- a/tests/test_linalg.py +++ b/tests/test_linalg.py @@ -386,6 +386,23 @@ def test_eig_arange(type, size): 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) + # and eigenvectors(v). + def assert_eigen_decomposition(self, a, w, v, rtol=1e-5, atol=1e-5): + a_ndim = a.ndim + if a_ndim == 2: + assert_allclose(a @ v, v @ inp.diag(w), rtol=rtol, atol=atol) + else: # a_ndim > 2 + if a_ndim > 3: + a = a.reshape(-1, *a.shape[-2:]) + w = w.reshape(-1, w.shape[-1]) + v = v.reshape(-1, *v.shape[-2:]) + for i in range(a.shape[0]): + assert_allclose( + a[i].dot(v[i]), w[i] * v[i], rtol=rtol, atol=atol + ) + @pytest.mark.parametrize( "func", [ @@ -413,11 +430,16 @@ def test_eigenvalues(self, func, shape, dtype, order): a_order = numpy.array(a, order=order) a_dp = inp.array(a, order=order) + # NumPy with OneMKL and with rocSOLVER sorts in ascending order, + # so w's should be directly comparable. + # 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, v = numpy.linalg.eigh(a_order) + w, _ = numpy.linalg.eigh(a_order) w_dp, v_dp = inp.linalg.eigh(a_dp) - assert_dtype_allclose(v_dp, v) + self.assert_eigen_decomposition(a_dp, w_dp, v_dp) else: # eighvalsh w = numpy.linalg.eigvalsh(a_order) From fa6078e6bad0c26a1a9e0c96b9f489b3f5e1a178 Mon Sep 17 00:00:00 2001 From: Vladislav Perevezentsev Date: Tue, 26 Mar 2024 19:16:29 +0100 Subject: [PATCH 26/26] Update test_eigenvalue in test_sycl_queue.py --- tests/test_sycl_queue.py | 14 +++++++++++++- 1 file changed, 13 insertions(+), 1 deletion(-) diff --git a/tests/test_sycl_queue.py b/tests/test_sycl_queue.py index c8dd69d594e5..dc531baa98bc 100644 --- a/tests/test_sycl_queue.py +++ b/tests/test_sycl_queue.py @@ -1264,7 +1264,19 @@ def test_eigenvalue(func, shape, device): dp_val, dp_vec = dpnp.linalg.eigh(dp_a) np_val, np_vec = numpy.linalg.eigh(a) - assert_allclose(dp_vec, np_vec, rtol=1e-05, atol=1e-05) + # Check the eigenvalue decomposition + if a.ndim == 2: + assert_allclose( + dp_a @ dp_vec, dp_vec @ dpnp.diag(dp_val), rtol=1e-5, atol=1e-5 + ) + else: # a.ndim == 3 + for i in range(a.shape[0]): + assert_allclose( + dp_a[i].dot(dp_vec[i]), + dp_val[i] * dp_vec[i], + rtol=1e-5, + atol=1e-5, + ) assert dp_vec.shape == np_vec.shape assert dp_vec.dtype == np_vec.dtype