///|
/// Add two matrices.  Invalid shapes produce a zero-sized matrix rather than
/// allowing a partially computed result into a filter.
pub fn Matrix::add(self : Matrix, other : Matrix) -> Matrix {
  if self.rows != other.rows || self.cols != other.cols {
    return Matrix::zeros(0, 0)
  }
  let result = Matrix::zeros(self.rows, self.cols)
  for i in 0.. Matrix {
  if self.rows != other.rows || self.cols != other.cols {
    return Matrix::zeros(0, 0)
  }
  let result = Matrix::zeros(self.rows, self.cols)
  for i in 0.. Matrix {
  let result = self.copy()
  for i in 0.. Matrix {
  if self.rows != other.rows || self.cols != other.cols {
    return Matrix::zeros(0, 0)
  }
  let result = Matrix::zeros(self.rows, self.cols)
  for i in 0.. Matrix {
  if self.cols != other.rows {
    return Matrix::zeros(0, 0)
  }
  let result = Matrix::zeros(self.rows, other.cols)
  for i in 0.. Array[Double] {
  if self.cols != vector.length() {
    return []
  }
  Array::makei(self.rows, i => {
    let mut total = 0.0
    for j in 0.. Matrix {
  let result = Matrix::zeros(self.cols, self.rows)
  for i in 0.. Double {
  if !self.is_square() {
    return 0.0
  }
  let mut result = 0.0
  for i in 0.. Matrix {
  let result = Matrix::zeros(left.length(), right.length())
  for i, a in left {
    for j, b in right {
      result.data[i * right.length() + j] = a * b
    }
  }
  result
}

///|
pub fn Matrix::symmetric_part(self : Matrix) -> Matrix {
  self.add(self.transpose()).scale(0.5)
}

///|
pub fn Matrix::skew_part(self : Matrix) -> Matrix {
  self.sub(self.transpose()).scale(0.5)
}

///|
pub fn Matrix::add_diagonal(self : Matrix, value : Double) -> Matrix {
  if !self.is_square() {
    return Matrix::zeros(0, 0)
  }
  let result = self.copy()
  for i in 0.. Matrix {
  if !self.is_square() || values.length() != self.rows {
    return Matrix::zeros(0, 0)
  }
  let result = self.copy()
  for i in 0.. Matrix {
  let rs = if row_start < 0 { 0 } else { row_start }
  let re = if row_end > self.rows { self.rows } else { row_end }
  let cs = if col_start < 0 { 0 } else { col_start }
  let ce = if col_end > self.cols { self.cols } else { col_end }
  if re <= rs || ce <= cs {
    return Matrix::zeros(0, 0)
  }
  let result = Matrix::zeros(re - rs, ce - cs)
  for i in 0.. Matrix {
  if row < 0 ||
    col < 0 ||
    row + source.rows > self.rows ||
    col + source.cols > self.cols {
    return Matrix::zeros(0, 0)
  }
  let result = self.copy()
  for i in 0.. Double {
  let mut total = 0.0
  for value in self.data {
    total = total + value * value
  }
  total.sqrt()
}

///|
/// Infinity norm: maximum absolute row sum.
pub fn Matrix::infinity_norm(self : Matrix) -> Double {
  let mut result = 0.0
  for i in 0.. result {
      result = row_sum
    }
  }
  result
}

///|
/// One norm: maximum absolute column sum.
pub fn Matrix::one_norm(self : Matrix) -> Double {
  let mut result = 0.0
  for j in 0.. result {
      result = column_sum
    }
  }
  result
}

///|
/// Construct a block diagonal matrix from independent square blocks.
pub fn Matrix::block_diagonal(blocks : Array[Matrix]) -> Matrix {
  let mut total = 0
  for block in blocks {
    if block.is_square() {
      total = total + block.rows
    }
  }
  let result = Matrix::zeros(total, total)
  let mut offset = 0
  for block in blocks {
    if block.is_square() {
      for i in 0..