ENH: linalg: use _umath_linalg for eigvalsh()

1 files changed, 20 insertions, 44 deletions

diff --git a/numpy/linalg/linalg.py b/numpy/linalg/linalg.py
index ae62d63fc..d36df8157 100644
--- a/numpy/linalg/linalg.py
+++ b/numpy/linalg/linalg.py
@@ -888,16 +888,16 @@ def eigvalsh(a, UPLO='L'):
 
     Parameters
     ----------
-    a : (M, M) array_like
+    a : (..., M, M) array_like
         A complex- or real-valued matrix whose eigenvalues are to be
         computed.
     UPLO : {'L', 'U'}, optional
-        Specifies whether the calculation is done with the lower triangular
-        part of `a` ('L', default) or the upper triangular part ('U').
+        Same as `lower`, wth 'L' for lower and 'U' for upper triangular.
+        Deprecated.
 
     Returns
     -------
-    w : (M,) ndarray
+    w : (..., M,) ndarray
         The eigenvalues, not necessarily ordered, each repeated according to
         its multiplicity.
 
@@ -915,9 +915,10 @@ def eigvalsh(a, UPLO='L'):
 
     Notes
     -----
-    This is a simple interface to the LAPACK routines dsyevd and zheevd
-    that sets those routines' flags to return only the eigenvalues of
-    real symmetric and complex Hermitian arrays, respectively.
+    Broadcasting rules apply, see the `numpy.linalg` documentation for
+    details.
+
+    The eigenvalues are computed using LAPACK routines _ssyevd, _heevd
 
     Examples
     --------
@@ -927,45 +928,20 @@ def eigvalsh(a, UPLO='L'):
     array([ 0.17157288+0.j,  5.82842712+0.j])
 
     """
-    UPLO = asbytes(UPLO)
+
+    extobj = get_linalg_error_extobj(
+        _raise_linalgerror_eigenvalues_nonconvergence)
+    if UPLO == 'L':
+        gufunc = _umath_linalg.eigvalsh_lo
+    else:
+        gufunc = _umath_linalg.eigvalsh_up
+
     a, wrap = _makearray(a)
-    _assertRank2(a)
-    _assertSquareness(a)
+    _assertNonEmpty(a)
+    _assertRankAtLeast2(a)
+    _assertNdSquareness(a)
     t, result_t = _commonType(a)
-    real_t = _linalgRealType(t)
-    a = _fastCopyAndTranspose(t, a)
-    a = _to_native_byte_order(a)
-    n = a.shape[0]
-    liwork = 5*n+3
-    iwork = zeros((liwork,), fortran_int)
-    if isComplexType(t):
-        lapack_routine = lapack_lite.zheevd
-        w = zeros((n,), real_t)
-        lwork = 1
-        work = zeros((lwork,), t)
-        lrwork = 1
-        rwork = zeros((lrwork,), real_t)
-        results = lapack_routine(_N, UPLO, n, a, n, w, work, -1,
-                                 rwork, -1, iwork, liwork,  0)
-        lwork = int(abs(work[0]))
-        work = zeros((lwork,), t)
-        lrwork = int(rwork[0])
-        rwork = zeros((lrwork,), real_t)
-        results = lapack_routine(_N, UPLO, n, a, n, w, work, lwork,
-                                rwork, lrwork, iwork, liwork,  0)
-    else:
-        lapack_routine = lapack_lite.dsyevd
-        w = zeros((n,), t)
-        lwork = 1
-        work = zeros((lwork,), t)
-        results = lapack_routine(_N, UPLO, n, a, n, w, work, -1,
-                                 iwork, liwork, 0)
-        lwork = int(work[0])
-        work = zeros((lwork,), t)
-        results = lapack_routine(_N, UPLO, n, a, n, w, work, lwork,
-                                 iwork, liwork, 0)
-    if results['info'] > 0:
-        raise LinAlgError('Eigenvalues did not converge')
+    w = gufunc(a.astype(t), extobj=extobj)
     return w.astype(result_t)
 
 def _convertarray(a):