///|
/// A dense matrix of `Double` in row-major order.
///
/// Element `(i, j)` lives at `data[i * cols + j]`. The arithmetic methods
/// return new matrices; `set` is the only writer.
#derive(Debug)
pub(all) struct Mat {
  rows : Int32
  cols : Int32
  data : FixedArray[Double]
}

///|
/// An all-zero `rows x cols` matrix.
pub fn Mat::zeros(rows : Int32, cols : Int32) -> Mat {
  { rows, cols, data: FixedArray::make(rows * cols, 0.0), }
}

///|
/// A `rows x cols` matrix with every entry `value`.
pub fn Mat::filled(rows : Int32, cols : Int32, value : Double) -> Mat {
  { rows, cols, data: FixedArray::make(rows * cols, value), }
}

///|
/// The `n x n` identity.
pub fn Mat::identity(n : Int32) -> Mat {
  let data = FixedArray::make(n * n, 0.0)
  for i in 0.. Mat {
  let n = values.len()
  let data = FixedArray::make(n * n, 0.0)
  for i in 0.. Mat {
  let rows = values.length()
  if rows == 0 {
    return Mat::zeros(0, 0)
  }
  let cols = values[0].length()
  let data = FixedArray::make(rows * cols, 0.0)
  for i in 0.. Mat {
  if data.length() != rows * cols {
    abort("Mat::from_fixed: storage does not match the shape")
  }
  { rows, cols, data, }
}

///|
/// The entry at `(i, j)`. An index outside the matrix aborts.
pub fn Mat::get(self : Mat, i : Int32, j : Int32) -> Double {
  self.data[i * self.cols + j]
}

///|
/// The entry at `(i, j)`, or `None` when the index is outside the matrix.
pub fn Mat::at(self : Mat, i : Int32, j : Int32) -> Double? {
  if i >= 0 && i < self.rows && j >= 0 && j < self.cols {
    Some(self.data[i * self.cols + j])
  } else {
    None
  }
}

///|
/// Overwrite the entry at `(i, j)`.
pub fn Mat::set(self : Mat, i : Int32, j : Int32, value : Double) -> Unit {
  self.data[i * self.cols + j] = value
}

///|
/// Add `value` to the entry at `(i, j)`.
pub fn Mat::add_to(self : Mat, i : Int32, j : Int32, value : Double) -> Unit {
  self.data[i * self.cols + j] += value
}

///|
/// The `i`-th row, copied out.
pub fn Mat::row(self : Mat, i : Int32) -> Vec {
  let data = FixedArray::make(self.cols, 0.0)
  for j in 0.. Vec {
  let data = FixedArray::make(self.rows, 0.0)
  for i in 0.. Unit {
  if values.len() != self.cols {
    abort("Mat::set_row: length mismatch")
  }
  for j in 0.. Unit {
  if values.len() != self.rows {
    abort("Mat::set_col: length mismatch")
  }
  for i in 0.. Mat {
  if row < 0 ||
    col < 0 ||
    height < 0 ||
    width < 0 ||
    row + height > self.rows ||
    col + width > self.cols {
    abort("Mat::block: block outside the matrix")
  }
  let out = Mat::zeros(height, width)
  for i in 0.. Mat {
  let data = FixedArray::make(self.rows * self.cols, 0.0)
  for i in 0.. Mat {
  let data = FixedArray::make(self.data.length(), 0.0)
  for k in 0.. Mat {
  if self.rows != other.rows || self.cols != other.cols {
    abort("Mat::add: shape mismatch")
  }
  let data = FixedArray::make(self.data.length(), 0.0)
  for k in 0.. Mat {
  if self.rows != other.rows || self.cols != other.cols {
    abort("Mat::sub: shape mismatch")
  }
  let data = FixedArray::make(self.data.length(), 0.0)
  for k in 0.. Mat {
  self.scale(-1.0)
}

///|
/// The product `self * other`.
pub fn Mat::mul(self : Mat, other : Mat) -> Mat {
  if self.cols != other.rows {
    abort("Mat::mul: shape mismatch")
  }
  let data = FixedArray::make(self.rows * other.cols, 0.0)
  for i in 0.. Vec {
  if v.len() != self.cols {
    abort("Mat::apply: length mismatch")
  }
  let data = FixedArray::make(self.rows, 0.0)
  for i in 0.. Vec {
  if v.len() != self.rows {
    abort("Mat::apply_transposed: length mismatch")
  }
  let data = FixedArray::make(self.cols, 0.0)
  for j in 0.. Double {
  let n = if self.rows < self.cols { self.rows } else { self.cols }
  let mut total = 0.0
  for i in 0.. Mat {
  if self.rows != right.rows {
    abort("Mat::augment: row count mismatch")
  }
  let out = Mat::zeros(self.rows, self.cols + right.cols)
  for i in 0.. Mat {
  if self.cols != below.cols {
    abort("Mat::stack: column count mismatch")
  }
  let out = Mat::zeros(self.rows + below.rows, self.cols)
  for i in 0.. Bool {
  for k in 0.. Bool {
  if self.rows != other.rows || self.cols != other.cols {
    return false
  }
  for k in 0.. tolerance {
      return false
    }
  }
  true
}

///|
/// The rows copied into growable arrays, for output and test expectations.
pub fn Mat::to_array_rows(self : Mat) -> Array[Array[Double]] {
  let out = Array::make(self.rows, Array::make(0, 0.0))
  for i in 0.. String {
  let lines = Array::make(self.rows, "")
  for i in 0..