///|
/// An LU factorisation with partial pivoting, `perm * A = lower * upper`, where
/// `lower` is unit lower triangular and `upper` is upper triangular.
///
/// `sign` is the parity of the row swaps, so `det = sign * product(upper)`. The
/// factorisation is kept so that several right-hand sides can reuse it, which
/// is what the estimators do inside a filter loop.
#derive(Debug)
pub(all) struct Lu {
  n : Int32
  lower : Mat
  upper : Mat
  perm : Array[Int32]
  sign : Double
}

///|
/// Factorise `a`, a square matrix, into `perm * a = lower * upper`.
///
/// Raises `Refused::Singular` at the first column whose pivot is zero;
/// `Refused::Dimension` when `a` is not square.
pub fn lu(a : Mat) -> Lu 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 work = Mat::zeros(n, n)
  for i in 0.. best {
        best = v
        pivot_row = i
      }
    }
    if best == 0.0 {
      raise Refused::Singular(at=k)
    }
    if pivot_row != k {
      for j in 0.. Mat raise Refused {
  if b.rows != self.n {
    raise Refused::Dimension(got=b.rows, want=self.n)
  }
  let width = b.cols
  let y = Mat::zeros(self.n, width)
  for i in 0..=..0 {
    for j in 0.. Double {
  let mut product = self.sign
  for i in 0..