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<sycl::event, sycl::event>
           const std::int8_t upper_lower,
           dpctl::tensor::usm_ndarray eig_vecs,
           dpctl::tensor::usm_ndarray eig_vals,
-          const std::vector<sycl::event> &depends);
+          const std::vector<sycl::event> &depends = {});
 
 extern void init_heevd_dispatch_table(void);
 } // namespace lapack
diff --git a/dpnp/linalg/dpnp_iface_linalg.py b/dpnp/linalg/dpnp_iface_linalg.py
index 07b2d46e0f7c..44301560c6bb 100644
--- a/dpnp/linalg/dpnp_iface_linalg.py
+++ b/dpnp/linalg/dpnp_iface_linalg.py
@@ -69,6 +69,7 @@
     "eig",
     "eigh",
     "eigvals",
+    "eigvalsh",
     "inv",
     "matrix_power",
     "matrix_rank",
@@ -247,6 +248,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
@@ -256,15 +270,13 @@ 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.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
     --------
@@ -282,20 +294,13 @@ def eigh(a, UPLO="L"):
     """
 
     dpnp.check_supported_arrays_type(a)
+    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'")
 
-    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)
 
 
@@ -327,6 +332,64 @@ def eigvals(input):
     return call_origin(numpy.linalg.eigvals, input)
 
 
+def eigvalsh(a, UPLO="L"):
+    """
+    eigvalsh(a, UPLO="L")
+
+    Compute 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, 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` : Return the eigenvalues and eigenvectors of a complex Hermitian
+                              (conjugate symmetric) or a real symmetric matrix.
+    :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
+    --------
+    >>> 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)
+
+    UPLO = UPLO.upper()
+    if UPLO not in ("L", "U"):
+        raise ValueError("UPLO argument must be 'L' or 'U'")
+
+    return dpnp_eigh(a, UPLO=UPLO, eigen_mode="N")
+
+
 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 3a1b565f02eb..2a50b67043d5 100644
--- a/dpnp/linalg/dpnp_utils_linalg.py
+++ b/dpnp/linalg/dpnp_utils_linalg.py
@@ -888,12 +888,14 @@ 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.
+    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
     for LAPACK library of OneMKL. Depending on input type of `a` array,
@@ -901,59 +903,76 @@ def dpnp_eigh(a, UPLO):
 
     """
 
-    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)
+    # get resulting type of arrays with eigenvalues and eigenvectors
+    v_type = _common_type(a)
+    w_type = _real_type(v_type)
 
-    # 'V' means both eigenvectors and eigenvalues will be calculated
-    jobz = _jobz["V"]
+    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
+
+    # `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
-    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
+    # 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"
+    )
+
+    a_sycl_queue = a.sycl_queue
+    a_order = "C" if a.flags.c_contiguous else "F"
 
     if a.ndim > 2:
-        w = dpnp.empty(
-            a.shape[:-1],
+        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])
+        a_usm_arr = dpnp.get_usm_ndarray(a)
+
+        # allocate a memory for dpnp array of eigenvalues
+        w = dpnp.empty_like(
+            a,
+            shape=orig_shape[:-1],
             dtype=w_type,
-            usm_type=a_usm_type,
-            sycl_queue=a_sycl_queue,
         )
+        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
-        op_count = a.shape[0]
-        if op_count == 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):
+        batch_size = a.shape[0]
+        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,
             )
 
+            # 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_lapack_ev[i], _ = getattr(li, lapack_func)(
+            ht_list_ev[2 * i + 1], _ = getattr(li, lapack_func)(
                 a_sycl_queue,
                 jobz,
                 uplo,
@@ -962,29 +981,43 @@ def dpnp_eigh(a, UPLO):
                 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)
+
+        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).reshape(orig_shape)
+            return w, v
+        return w
 
-        # combine the list of eigenvectors into a single array
-        v = dpnp.array(eig_vecs, order=a_order)
-        return w, v
     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)
+        a_usm_arr = dpnp.get_usm_ndarray(a)
+        ht_list_ev = []
+        copy_ev = dpctl.SyclEvent()
 
-        # 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 can 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(
-            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
@@ -996,8 +1029,9 @@ def dpnp_eigh(a, UPLO):
             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(
@@ -1006,14 +1040,13 @@ def dpnp_eigh(a, UPLO):
                 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_inv_batched(a, res_type):
diff --git a/tests/helper.py b/tests/helper.py
index 2a2873afdcea..ba6e32a247c8 100644
--- a/tests/helper.py
+++ b/tests/helper.py
@@ -166,6 +166,62 @@ def get_all_dtypes(
     return dtypes
 
 
+def generate_random_numpy_array(
+    shape, dtype=None, hermitian=False, seed_value=None
+):
+    """
+    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
+    seed_value : int, optional
+        The seed value to initialize the random number generator.
+        Default : None
+
+    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, the Hermitian or
+    symmetric property is ensured for each 2D sub-array.
+
+    """
+
+    numpy.random.seed(seed_value)
+
+    a = numpy.random.randn(*shape).astype(dtype)
+    if numpy.issubdtype(a.dtype, numpy.complexfloating):
+        numpy.random.seed(seed_value)
+        a += 1j * numpy.random.randn(*shape)
+
+    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
+
+
 def is_cpu_device(device=None):
     """
     Return True if a test is running on CPU device, False otherwise.
diff --git a/tests/skipped_tests_gpu_no_fp64.tbl b/tests/skipped_tests_gpu_no_fp64.tbl
index fc591b36bd60..13eaf9f39211 100644
--- a/tests/skipped_tests_gpu_no_fp64.tbl
+++ b/tests/skipped_tests_gpu_no_fp64.tbl
@@ -20,9 +20,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/test_linalg.py b/tests/test_linalg.py
index 6b9800a8d9b3..53668819a10a 100644
--- a/tests/test_linalg.py
+++ b/tests/test_linalg.py
@@ -15,6 +15,7 @@
 
 from .helper import (
     assert_dtype_allclose,
+    generate_random_numpy_array,
     get_all_dtypes,
     get_complex_dtypes,
     get_float_complex_dtypes,
@@ -383,47 +384,93 @@ 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.SyclDevice().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 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",
+        [
+            "eigh",
+            "eigvalsh",
+        ],
     )
-    symm = symm_orig
-    dpnp_symm_orig = inp.array(symm)
-    dpnp_symm = dpnp_symm_orig
+    @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)
 
-    dpnp_val, dpnp_vec = inp.linalg.eigh(dpnp_symm)
-    np_val, np_vec = numpy.linalg.eigh(symm)
+        # 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, _ = numpy.linalg.eigh(a_order)
+            w_dp, v_dp = inp.linalg.eigh(a_dp)
 
-    # DPNP sort val/vec by abs value
-    vvsort(dpnp_val, dpnp_vec, size, inp)
+            self.assert_eigen_decomposition(a_dp, w_dp, v_dp)
 
-    # NP sort val/vec by abs value
-    vvsort(np_val, np_vec, size, numpy)
+        else:  # eighvalsh
+            w = numpy.linalg.eigvalsh(a_order)
+            w_dp = inp.linalg.eigvalsh(a_dp)
 
-    # 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]
+        assert_dtype_allclose(w_dp, w)
 
-    assert_array_equal(symm_orig, symm)
-    assert_array_equal(dpnp_symm_orig, dpnp_symm)
+    @pytest.mark.parametrize(
+        "func",
+        [
+            "eigh",
+            "eigvalsh",
+        ],
+    )
+    def test_eigenvalue_errors(self, func):
+        a_dp = inp.array([[1, 3], [3, 2]], dtype="float32")
 
-    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
+        # unsupported type
+        a_np = inp.asnumpy(a_dp)
+        dpnp_func = getattr(inp.linalg, func)
+        assert_raises(TypeError, dpnp_func, a_np)
 
-    assert_allclose(dpnp_val, np_val, rtol=1e-05, atol=1e-04)
-    assert_allclose(dpnp_vec, np_vec, rtol=1e-05, atol=1e-04)
+        # a.ndim < 2
+        a_dp_ndim_1 = a_dp.flatten()
+        assert_raises(inp.linalg.LinAlgError, dpnp_func, a_dp_ndim_1)
+
+        # a is not square
+        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")
 
 
 @pytest.mark.parametrize("type", get_all_dtypes(no_bool=True, no_complex=True))
@@ -1127,11 +1174,6 @@ def test_norm_error(self):
 
 
 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
@@ -1156,9 +1198,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 seed_value=81 to prevent
+        # random generation of the input singular matrix
+        a = generate_random_numpy_array(shape, dtype, seed_value=81)
         ia = inp.array(a)
 
         if mode == "r":
@@ -1228,7 +1270,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
@@ -1484,11 +1526,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):
@@ -1586,11 +1623,11 @@ 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
-
+        # 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 = inp.array(a)
 
         if compute_vt:
@@ -1628,11 +1665,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):
@@ -1668,9 +1700,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 seed_value=81 to prevent
+        # random generation of the input singular matrix
+        a = generate_random_numpy_array(shape, dtype, seed_value=81)
         a_dp = inp.array(a)
 
         B = numpy.linalg.pinv(a)
@@ -1695,11 +1727,11 @@ 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
-
+        # 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_dp = inp.array(a)
 
         B = numpy.linalg.pinv(a, hermitian=True)
@@ -1735,7 +1767,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 1043d4444e93..f080f18590b6 100644
--- a/tests/test_sycl_queue.py
+++ b/tests/test_sycl_queue.py
@@ -11,7 +11,12 @@
 from dpnp.dpnp_array import dpnp_array
 from dpnp.dpnp_utils import get_usm_allocations
 
-from .helper import assert_dtype_allclose, get_all_dtypes, is_win_platform
+from .helper import (
+    assert_dtype_allclose,
+    generate_random_numpy_array,
+    get_all_dtypes,
+    is_win_platform,
+)
 
 list_of_backend_str = [
     "host",
@@ -1217,43 +1222,79 @@ 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(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
-
-    dpnp_val, dpnp_vec = dpnp.linalg.eigh(dpnp_data)
-    numpy_val, numpy_vec = numpy.linalg.eigh(numpy_data)
-
-    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
-
+    # 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
+
+    if func == "eigh":
+        dp_val, dp_vec = dpnp.linalg.eigh(dp_a)
+        np_val, np_vec = numpy.linalg.eigh(a)
+
+        # 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
+
+        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
+    assert dp_val.dtype == np_val.dtype
+
+    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(
@@ -1944,13 +1985,12 @@ 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)
-    if hermitian:
-        a_np = numpy.random.randn(*shape) + 1j * numpy.random.randn(*shape)
-        a_np = numpy.conj(a_np.T) @ a_np
-    else:
-        a_np = numpy.random.randn(*shape)
-
+    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, seed_value=81
+    )
     a_dp = dpnp.array(a_np, device=device)
 
     if rcond_as_array:
diff --git a/tests/test_usm_type.py b/tests/test_usm_type.py
index 2eccb30ba283..83cf35198ab8 100644
--- a/tests/test_usm_type.py
+++ b/tests/test_usm_type.py
@@ -9,7 +9,7 @@
 import dpnp as dp
 from dpnp.dpnp_utils import get_usm_allocations
 
-from .helper import assert_dtype_allclose
+from .helper import assert_dtype_allclose, generate_random_numpy_array
 
 list_of_usm_types = ["device", "shared", "host"]
 
@@ -810,6 +810,45 @@ def test_where(usm_type):
     assert result.usm_type == 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(func, shape, usm_type):
+    a_np = generate_random_numpy_array(shape, hermitian=True)
+    a = dp.array(a_np, 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
 )
@@ -1047,14 +1086,9 @@ def test_matrix_rank(data, tol, usm_type):
     ],
 )
 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
diff --git a/tests/third_party/cupy/linalg_tests/test_eigenvalue.py b/tests/third_party/cupy/linalg_tests/test_eigenvalue.py
index b620bd39e984..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,13 +20,12 @@ def _get_hermitian(xp, a, UPLO):
         }
     )
 )
-class TestEigenvalue(unittest.TestCase):
+class TestEigenvalue:
     @testing.for_all_dtypes()
     @testing.numpy_cupy_allclose(
         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 +62,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 +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, contiguous_check=False)
+    @testing.numpy_cupy_allclose(rtol=1e-3, atol=1e-4)
     def test_eigh_complex_batched(self, xp, dtype):
         a = xp.array(
             [
@@ -113,9 +112,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 +123,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,7 +140,6 @@ 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)
     def test_eigvalsh_complex(self, xp, dtype):
@@ -149,7 +149,6 @@ 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)
     def test_eigvalsh_complex_batched(self, xp, dtype):
@@ -171,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):
@@ -179,11 +178,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)
 
@@ -196,16 +194,15 @@ 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)
-            with pytest.raises((numpy.linalg.LinAlgError, ValueError)):
+            with pytest.raises(xp.linalg.LinAlgError):
                 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)