///|
/// Interior row-major COO, with row/column indices in 0.. InteriorMatrix raise SolverError {
  let g = self.grid
  if max_entries < 0 || g.unknowns() > max_entries / 5 {
    raise Budget("COO export exceeds entry budget")
  }
  let entries : Array[(Int32, Int32, Double)] = []
  for j in 1.. 1 {
        entries.push((row, row - 1, -self.west[p]))
      }
      if i < g.nx - 1 {
        entries.push((row, row + 1, -self.east[p]))
      }
      if j > 1 {
        entries.push((row, row - (g.nx - 1), -self.south[p]))
      }
      if j < g.ny - 1 {
        entries.push((row, row + (g.nx - 1), -self.north[p]))
      }
    }
  }
  { grid: g, entries, }
}

///|
pub fn InteriorMatrix::entries(
  self : InteriorMatrix,
) -> Array[(Int32, Int32, Double)] {
  self.entries.copy()
}

///|
pub fn InteriorMatrix::dimension(self : InteriorMatrix) -> Int32 {
  self.grid.unknowns()
}

///|
pub fn Field::interior_values(self : Field) -> Array[Double] {
  let g = self.grid
  let result = []
  for j in 1.. Field raise SolverError {
  if values.length() != grid.unknowns() {
    raise Shape("interior vector length mismatch")
  }
  let result = Array::make(grid.size(), 0.0)
  for j in 1.. Array[Double] raise SolverError {
  if values.length() != self.dimension() {
    raise Shape("matrix vector length mismatch")
  }
  for v in values {
    require_finite(v, "matrix input")
  }
  let result = Array::make(self.dimension(), 0.0)
  for entry in self.entries {
    let (row, column, value) = entry
    result[row] = result[row] + value * values[column]
  }
  for v in result {
    require_finite(v, "matrix result")
  }
  result
}