///|
pub(all) struct RelaxationResult {
  values : Array[Double]
  residuals : Array[Double]
  iterations : Int
  converged : Bool
} derive(Debug, ToJson)

///|
pub fn vector_distance(a : ArrayView[Double], b : ArrayView[Double]) -> Double {
  let count = if a.length() < b.length() { a.length() } else { b.length() }
  let mut total = 0.0
  for i in 0.. RelaxationResult {
  let values = initial.copy()
  let residuals : Array[Double] = []
  let mut residual = vector_distance(values, target)
  let mut iterations = 0
  while iterations < max_iterations && residual > 1.0e-12 {
    residuals.push(residual)
    for i in 0.. Double {
  l2_norm(residual)
}

///|
pub fn relative_residual(
  values : ArrayView[Double],
  target : ArrayView[Double],
) -> Double {
  safe_ratio(vector_distance(values, target), l2_norm(target), 0.0)
}

///|
pub fn diagonal_solve(
  diagonal : ArrayView[Double],
  right_hand_side : ArrayView[Double],
) -> Array[Double] {
  let output = zeros(diagonal.length())
  for i in 0.. Array[Double] {
  let output = zeros(rhs.length())
  for i in 0.. 0 && i - 1 < lower.length() {
      value = value + lower[i - 1] * solution[i - 1]
    }
    if i < upper.length() && i + 1 < solution.length() {
      value = value + upper[i] * solution[i + 1]
    }
    output[i] = value - rhs[i]
  }
  output
}

///|
pub fn tridiagonal_jacobi(
  diagonal : ArrayView[Double],
  lower : ArrayView[Double],
  upper : ArrayView[Double],
  rhs : ArrayView[Double],
  iterations : Int,
) -> RelaxationResult {
  let solution = zeros(rhs.length())
  let residuals : Array[Double] = []
  let mut step = 0
  while step < iterations {
    let residual = tridiagonal_residual(diagonal, lower, upper, solution, rhs)
    residuals.push(residual_norm(residual))
    for i in 0.. 0 && i - 1 < lower.length() {
        off = off + lower[i - 1] * solution[i - 1]
      }
      if i < upper.length() && i + 1 < solution.length() {
        off = off + upper[i] * solution[i + 1]
      }
      solution[i] = safe_ratio(rhs[i] - off, diagonal[i], solution[i])
    }
    step = step + 1
  }
  {
    values: solution,
    residuals,
    iterations: step,
    converged: residuals.length() > 0 &&
    residuals[residuals.length() - 1] < 1.0e-8,
  }
}

///|
pub fn solver_convergence_rate(residuals : ArrayView[Double]) -> Double {
  if residuals.length() < 2 {
    0.0
  } else {
    safe_ratio(residuals[residuals.length() - 1], residuals[0], 0.0)
  }
}

///|
pub fn solver_is_decreasing(residuals : ArrayView[Double]) -> Bool {
  let mut result = true
  for i in 1.. residuals[i - 1] {
      result = false
    }
  }
  result
}

///|
pub fn solver_residual_csv(residuals : ArrayView[Double]) -> String {
  let output = StringBuilder()
  output.write_string("iteration,residual\n")
  for i in 0..