///|
/// Determinant computed with partial pivoting.  A zero result means that the
/// matrix is singular or not square.
pub fn Matrix::determinant(self : Matrix) -> Double {
  if !self.is_square() {
    return 0.0
  }
  if self.rows == 0 {
    return 1.0
  }
  let work = self.copy()
  let mut sign = 1.0
  let mut determinant = 1.0
  for column in 0.. magnitude {
        magnitude = candidate
        pivot = row
      }
    }
    if magnitude < 0.000000000001 {
      return 0.0
    }
    if pivot != column {
      for j in column.. ignore
        work.set(pivot, j, temp) |> ignore
      }
      sign = -sign
    }
    let diagonal = work.get(column, column)
    determinant = determinant * diagonal
    for row in (column + 1).. ignore
      }
    }
  }
  determinant * sign
}

///|
/// Solve `A x = b` with partial pivoting.
pub fn Matrix::solve(self : Matrix, rhs : Array[Double]) -> MatrixSolveResult {
  if !self.is_square() || rhs.length() != self.rows {
    return InvalidShape
  }
  let n = self.rows
  let augmented = Matrix::zeros(n, n + 1)
  for i in 0.. magnitude {
        magnitude = candidate
        pivot = row
      }
    }
    if magnitude < 0.000000000001 {
      return Singular
    }
    if pivot != column {
      for j in column.. ignore
        augmented.set(pivot, j, temp) |> ignore
      }
    }
    let diagonal = augmented.get(column, column)
    for j in column.. ignore
    }
    for row in 0.. ignore
          }
        }
      }
    }
  }
  Solved(Array::makei(n, i => augmented.get(i, n)))
}

///|
/// Invert a matrix with the same pivoting strategy used by `solve`.
pub fn Matrix::inverse(self : Matrix) -> Matrix? {
  if !self.is_square() {
    return None
  }
  let n = self.rows
  let augmented = Matrix::zeros(n, n * 2)
  for i in 0.. magnitude {
        magnitude = candidate
        pivot = row
      }
    }
    if magnitude < 0.000000000001 {
      return None
    }
    if pivot != column {
      for j in 0.. ignore
        augmented.set(pivot, j, temp) |> ignore
      }
    }
    let diagonal = augmented.get(column, column)
    for j in 0.. ignore
    }
    for row in 0.. ignore
          }
        }
      }
    }
  }
  let result = Matrix::zeros(n, n)
  for i in 0.. Matrix? {
  if !self.is_square() {
    return None
  }
  let n = self.rows
  let result = Matrix::zeros(n, n)
  for i in 0.. ignore
      } else {
        let diagonal = result.get(j, j)
        if diagonal.abs() < 0.000000000001 {
          return None
        }
        result.set(i, j, sum / diagonal) |> ignore
      }
    }
  }
  Some(result)
}

///|
/// Add a small diagonal jitter until Cholesky succeeds or the budget is
/// exhausted.  This is useful for covariance matrices assembled from noisy
/// samples.
pub fn Matrix::regularized_cholesky(
  self : Matrix,
  initial_jitter : Double,
  attempts : Int,
) -> (Matrix, Double)? {
  if !self.is_square() {
    return None
  }
  let mut jitter : Double = if initial_jitter <= 0.0 {
    0.000000001
  } else {
    initial_jitter
  }
  let tries = if attempts < 1 { 1 } else { attempts }
  for _ in 0.. return Some((factor, jitter))
      None => jitter = jitter * 10.0
    }
  }
  None
}

///|
/// Rank estimate based on pivot magnitudes.
pub fn Matrix::rank(self : Matrix, tolerance : Double) -> Int {
  if self.is_empty() {
    return 0
  }
  let work = self.copy()
  let threshold : Double = if tolerance <= 0.0 {
    0.0000000001
  } else {
    tolerance
  }
  let mut row = 0
  let mut rank = 0
  for column in 0..= self.rows {
      break
    }
    let mut pivot = row
    let mut magnitude = work.get(row, column).abs()
    for candidate in (row + 1).. magnitude {
        magnitude = current
        pivot = candidate
      }
    }
    if magnitude <= threshold {
      continue
    }
    if pivot != row {
      for j in column.. ignore
        work.set(pivot, j, temp) |> ignore
      }
    }
    for candidate in (row + 1).. ignore
      }
    }
    rank = rank + 1
    row = row + 1
  }
  rank
}

///|
/// Jacobi sweeps for a symmetric matrix.  It is intentionally bounded: a
/// diagnostic should always return even for malformed field data.
pub fn Matrix::jacobi_eigenvalues(self : Matrix, sweeps : Int) -> Array[Double] {
  if !self.is_square() {
    return []
  }
  let work = self.symmetric_part()
  let count = if sweeps < 1 { 1 } else { sweeps }
  for _ in 0.. largest {
          largest = magnitude
          p = i
          q = j
        }
      }
    }
    if largest < 0.000000000001 {
      break
    }
    let app = work.get(p, p)
    let aqq = work.get(q, q)
    let apq = work.get(p, q)
    let angle = 0.5 * (aqq - app) / apq
    let sign = if angle >= 0.0 { 1.0 } else { -1.0 }
    let tangent = sign / (angle.abs() + (1.0 + angle * angle).sqrt())
    let cosine = 1.0 / (1.0 + tangent * tangent).sqrt()
    let sine = tangent * cosine
    for k in 0.. ignore
      work.set(k, q, sine * kip + cosine * kiq) |> ignore
    }
    for k in 0.. ignore
      work.set(q, k, sine * pk + cosine * qk) |> ignore
    }
  }
  Array::makei(work.rows, i => work.get(i, i))
}