// matdet  matrix determinant
func matdet(_ a: [[Double]]) -> Double {
  let n=a.count
  var D = 1.0                 // determinant
  var B = [[Double]](repeating:[Double](repeating:0.0,count:n),count:n)
                                        // working matrix
  var row = [Int](repeating:0,count:n)  // row interchange indicies
  var hold = 0
  var I_pivot = 0                       // pivot indicies
  var pivot = 0.0                       // pivot element value
  var abs_pivot = 0.0

  if a[0].count != n {
    print("Error in Matrix.determinant, inconsistent array sizes.")
  }
  // build working matrix
  for i in 0..<n {
    for j in 0..<n {
      B[i][j] = a[i][j]
    }
  }
  // set up row interchange vectors
  for k in 0..<n {
    row[k] = k
  }
  // begin main reduction loop
  for k in 0..<n-1 {
    // find largest element for pivot
    pivot = B[row[k]][k]
    abs_pivot = abs(pivot)
    I_pivot = k
    for i in k..<n {
      if abs(B[row[i]][k]) > abs_pivot {
        I_pivot = i
        pivot = B[row[i]][k]
        abs_pivot = abs(pivot)
      }
    }
    // have pivot, interchange row indicies
    if I_pivot != k {
      hold = row[k]
      row[k] = row[I_pivot]
      row[I_pivot] = hold
      D = -D
    }
    // check for near singular
    if abs_pivot < 1.0E-20 {
        return abs_pivot
    }
    else {
      D = D * pivot
      // reduce about pivot
      for j in k+1..<n {
          B[row[k]][j] = B[row[k]][j] / B[row[k]][k]
      }
      //  inner reduction loop
      for i in 0..<n {
        if i != k {
          for j in k+1..<n {
            B[row[i]][j] = B[row[i]][j] - B[row[i]][k] * B[row[k]][j]
          }
        }
      }
    }
    //  finished inner reduction
  }
  // end of main reduction loop
  return D * B[row[n-1]][n-1]
} // end matdet