///|
/// Descending eigenvalues with corresponding unit eigenvectors stored as rows.
pub struct Spectrum {
  values : Array[Double]
  vectors : Array[Array[Double]]
  residual : Double
}

///|
pub fn Spectrum::eigenvalues(self : Spectrum) -> Array[Double] {
  self.values.copy()
}

///|
pub fn Spectrum::eigenvectors(self : Spectrum) -> Array[Array[Double]] {
  self.vectors.map(r => r.copy())
}

///|
/// Largest absolute entry of A*v-lambda*v, in original matrix units.
pub fn Spectrum::max_residual(self : Spectrum) -> Double {
  self.residual
}

///|
/// Symmetric cyclic Jacobi, scaled before iteration; no silent nonconvergence.
pub fn symmetric_eigen(
  matrix : Array[Array[Double]],
) -> Spectrum raise CodaError {
  let n = matrix.length()
  if n < 1 || n > 63 {
    raise Budget("eigen dimension must be 1..63")
  }
  let mut scale = 0.0
  for row in matrix {
    if row.length() != n {
      raise Shape("eigen input must be square")
    }
    real_vector(row)
    for x in row {
      if x.abs() > scale {
        scale = x.abs()
      }
    }
  }
  if scale == 0.0 {
    scale = 1.0
  }
  for i in 0.. 1.0e-12 {
        raise Invalid("eigen input is not symmetric")
      }
    }
  }
  let a = Array::makei(n, i => {
    Array::makei(n, j => (matrix[i][j] / scale + matrix[j][i] / scale) * 0.5)
  })
  let v = Array::makei(n, i => {
    Array::makei(n, j => if i == j { 1.0 } else { 0.0 })
  })
  let mut converged = false
  for _sweep in 0..<100 {
    let mut largest = 0.0
    for p in 0.. largest {
          largest = apq.abs()
        }
        if apq.abs() > 1.0e-15 {
          let tau = (a[q][q] - a[p][p]) / (2.0 * apq)
          let sign = if tau < 0.0 { -1.0 } else { 1.0 }
          let t = sign / (tau.abs() + (1.0 + tau * tau).sqrt())
          let c = 1.0 / (1.0 + t * t).sqrt()
          let s = t * c
          let app = a[p][p]
          let aqq = a[q][q]
          a[p][p] = app - t * apq
          a[q][q] = aqq + t * apq
          a[p][q] = 0.0
          a[q][p] = 0.0
          for k in 0.. i)
  for i in 0.. a[order[best]][order[best]] {
        best = j
      }
    }
    let old = order[i]
    order[i] = order[best]
    order[best] = old
  }
  let values = order.map(i => checked(a[i][i] * scale))
  let vectors = order.map(i => {
    let row = Array::makei(n, j => v[j][i])
    let mut pivot = 0
    for j in 1.. row[pivot].abs() {
        pivot = j
      }
    }
    if row[pivot] < 0.0 {
      for j in 0.. residual {
        residual = r.abs()
      }
    }
  }
  if residual > 1.0e-9 {
    raise Numerical("eigen residual exceeds tolerance")
  }
  { values, vectors, residual: checked(residual * scale), }
}