docs: def eigh(a, eigvec=True, rcond=None):
      """ Eigenvalues and eigenvectors of symmetric matrix ``a``.
  
      Args:
          a: Two-dimensional, square Hermitian matrix/array of numbers
              and/or :class:`gvar.GVar`\s. Array elements must be
              real-valued if `gvar.GVar`\s are involved (i.e., symmetric
              matrix).
          eigvec (bool): If ``True``  (default), method returns a tuple
              of arrays ``(val, vec)`` where ``val[i]`` are the
              eigenvalues of ``a`` (in ascending order), and ``vec[:, i]``
              are the corresponding eigenvectors of ``a``. Only ``val`` is
              returned if ``eigvec=False``.
          rcond (float): Eigenvalues whose difference is smaller than
              ``rcond`` times their sum are assumed to be degenerate
              (and ignored) when computing variances for the eigvectors.
              Default (``rcond=None``) is ``max(M,N)`` times machine precision.
  
      Returns:
          Tuple ``(val,vec)`` of eigenvalues and eigenvectors of
          matrix ``a`` if parameter ``eigvec==True`` (default).
          The eigenvalues ``val[i]`` are in ascending order and
          ``vec[:, i]`` are the corresponding eigenvalues. Only
          the eigenvalues ``val`` are returned if ``eigvec=False``.
  
      Raises:
          ValueError: If matrix is not square and two-dimensional.
      """
      a = numpy.asarray(a)
      if a.dtype != object:
          val, vec = numpy.linalg.eigh(a)
          return (val, vec) if eigvec else val
      amean = gvar.mean(a)
      if amean.ndim != 2 or amean.shape[0] != amean.shape[1]:
          raise ValueError('bad matrix shape: ' + str(a.shape))
      if rcond is None:
          rcond = numpy.finfo(float).eps * max(a.shape)
      da = a - amean
      val0, vec0 = numpy.linalg.eigh(amean)
      val = val0 + [
          vec0[:, i].conjugate().dot(da.dot(vec0[:, i])) for i in range(vec0.shape[1])
          ]
      if eigvec == True:
          if vec0.dtype == complex:
              raise ValueError('cannot evaluate eigenvectors when a is complex')
          vec = numpy.array(vec0, dtype=object)
          for i in range(len(val)):
              for j in range(len(val)):
                  dval = val0[i] - val0[j]
                  if abs(dval) < rcond * abs(val0[j] + val0[i]) or dval == 0.0:
                      continue
                  vec[:, i] += vec0[:, j] * (
                      vec0[:, j].dot(da.dot(vec0[:, i])) / dval
                      )
          return val, vec
      else:
          return val