///|
fn matrix_shape(a : Array[Array[Double]]) -> (Int32, Int32) raise DesignError {
  if a.length() == 0 || a.length() > 4096 {
    raise Budget("matrix row count must be 1..4096")
  }
  let cols = a[0].length()
  if cols == 0 || cols > 128 {
    raise Budget("matrix column count must be 1..128")
  }
  for row in a {
    if row.length() != cols {
      raise Invalid("ragged matrix")
    }
    for x in row {
      require_number(x)
    }
  }
  (a.length(), cols)
}

///|
fn zeros(rows : Int32, cols : Int32) -> Array[Array[Double]] {
  [
    for _ in 0.. Array::make(cols, 0.0)
  ]
}

///|
fn tolerance(tol : Double) -> Unit raise DesignError {
  if !finite(tol) || tol < 1.0e-14 || tol > 0.01 {
    raise Invalid("relative tolerance must be 1e-14..0.01")
  }
}

///|
pub fn transpose(
  a : Array[Array[Double]],
) -> Array[Array[Double]] raise DesignError {
  let (n, p) = matrix_shape(a)
  let result = zeros(p, n)
  for i in 0.. Array[Array[Double]] raise DesignError {
  let (n, p) = matrix_shape(a)
  let (m, q) = matrix_shape(b)
  if p != m {
    raise Invalid("matrix multiplication shape mismatch")
  }
  if q > 128 {
    raise Budget("product exceeds column budget")
  }
  let result = zeros(n, q)
  for i in 0.. Int32 raise DesignError {
  tolerance(tol)
  let (n, p) = matrix_shape(a)
  let work = [ for row in a => row.copy() ]
  // Equilibrate columns: the threshold must not depend on a user's physical units.
  for j in 0..
 scale {
        scale = work[i][j].abs()
      }
    }
    if scale > 0.0 {
      for i in 0.. work[pivot][col].abs() {
        pivot = i
      }
    }
    if work[pivot][col].abs() <= tol {
      continue
    }
    let swap = work[rank]
    work[rank] = work[pivot]
    work[pivot] = swap
    for i in (rank + 1).. Array[Array[Double]] raise DesignError {
  let (n, p) = matrix_shape(a)
  if weights.length() != 0 && weights.length() != n {
    raise Invalid("weight count mismatch")
  }
  let result = zeros(p, p)
  for i in 0..