///|
/// The eigenvalues and eigenvectors of a symmetric matrix.
///
/// `values` is ascending and `vectors` has the matching eigenvector in its
/// matching column, so `a * v[:, i] = values[i] * v[:, i]`.
#derive(Debug)
pub(all) struct Eigh {
  values : Vec
  vectors : Mat
}

///|
/// The symmetric eigendecomposition of `a` by cyclic Jacobi rotations.
///
/// The rotations are stable for the covariance-shaped matrices this library
/// hands it, and they need no shift strategy. Only the symmetric part of `a` is
/// used, so a small skew from rounding does not matter. Raises
/// `Refused::Dimension` when `a` is not square and `Refused::Exceeded` when the
/// sweep limit is reached.
pub fn eigh(a : Mat, sweeps~ : Int32 = 64) -> Eigh raise Refused {
  let n = a.rows
  if a.rows != a.cols {
    raise Refused::Dimension(got=a.cols, want=a.rows)
  }
  if n == 0 {
    raise Refused::Empty
  }
  let m = Mat::zeros(n, n)
  for i in 0..= 0.0 {
          1.0 / (theta + (1.0 + theta * theta).sqrt())
        } else {
          -1.0 / (-theta + (1.0 + theta * theta).sqrt())
        }
        let c = 1.0 / (1.0 + t * t).sqrt()
        let s = t * c
        for k in 0..= sweeps {
    raise Refused::Exceeded(limit=sweeps, got=sweeps)
  }
  let values = Vec::zeros(n)
  for i in 0..