///|
pub fn inverse(
  a : Array[Array[Double]],
  tol~ : Double = 1.0e-10,
) -> Array[Array[Double]] raise DesignError {
  tolerance(tol)
  let (n, p) = matrix_shape(a)
  if n != p {
    raise Invalid("inverse requires square matrix")
  }
  let work = [ for row in a => row.copy() ]
  let result = zeros(n, n)
  for i in 0.. scale {
        scale = x.abs()
      }
    }
    if scale == 0.0 {
      raise Singular("zero matrix row")
    }
    for j in 0.. work[pivot][col].abs() {
        pivot = i
      }
    }
    if work[pivot][col].abs() <= tol {
      raise Singular("matrix pivot below tolerance")
    }
    let row = work[col]
    work[col] = work[pivot]
    work[pivot] = row
    let right = result[col]
    result[col] = result[pivot]
    result[pivot] = right
    let value = work[col][col]
    for j in 0.. Double? raise DesignError {
  tolerance(tol)
  let (n, p) = matrix_shape(a)
  if n != p {
    raise Invalid("positive determinant requires square matrix")
  }
  let scales = Array::make(n, 0.0)
  for i in 0.. tol * magnitude {
        raise Invalid("matrix must be symmetric")
      }
    }
  }
  let l = zeros(n, n)
  let mut value = 0.0
  for i in 0..