// matrix2.swift
// Numerical matrix examples
#if os(OSX) || os(iOS)
import Darwin
#elseif os(Linux)
import Glibc
#endif

// add two n by m matrices
func matadd(_ a: [[Double]], _ b: [[Double]]) -> [[Double]] {
  let nia = a.count
  let nja = a[0].count
  if nia != b.count || nja != b[0].count {
    return [[Double]](repeating:[Double](repeating:0.0,count:nja),count:nia)
  }
  var sum = [[Double]](repeating:[Double](repeating:0.0,count:nja),count:nia)
  for i in 0..<nia {
    for j in 0..<nja {
      sum[i][j] = a[i][j]+b[i][j]
    }
  }
  return sum
} // end matadd

// matinv  matrix invert
func matinv(_ a: [[Double]]) -> [[Double]] {
  let n = a.count
  var inv = [[Double]](repeating:[Double](repeating:0.0,count:n),count:n)
  for i in 0..<n {
    for j in 0..<n {
      inv[i][j] = a[i][j]
    }
  }

  var row = [Int](repeating:0,count:n)
  var col = [Int](repeating:0,count:n)
  var temp = [Double](repeating:0,count:n)
  var hold = 0
  var I_pivot = 0
  var J_pivot = 0
  var pivot = 0.0
  var abs_pivot = 0.0

  // set up row and column interchange vectors
  for k in 0..<n {
    row[k] = k
    col[k] = k
  }
  // begin main reduction loop
  for k in 0..<n {
    // find largest element for pivot
    pivot = inv[row[k]][col[k]]
    I_pivot = k
    J_pivot = k
    for i in k..<n {
      for j in k..<n {
        abs_pivot = abs(pivot)
        if abs(inv[row[i]][col[j]]) > abs_pivot {
          I_pivot = i
          J_pivot = j
          pivot = inv[row[i]][col[j]]
        }
      }
    }
    if abs(pivot) < 1.0E-10 {
      print("Matrix is singular !")
      return [[Double]](repeating:[Double](repeating:0.0,count:n),count:n)
    }
    hold = row[k]
    row[k] = row[I_pivot]
    row[I_pivot] = hold
    hold = col[k]
    col[k] = col[J_pivot]
    col[J_pivot] = hold
    // reduce about pivot
    inv[row[k]][col[k]] = 1.0 / pivot
    for j in 0..<n {
      if j != k {
        inv[row[k]][col[j]] = inv[row[k]][col[j]] * inv[row[k]][col[k]]
      }
    }
    // inner reduction loop
    for i in 0..<n {
      if k != i {
        for j in 0..<n {
          if k != j {
            inv[row[i]][col[j]] = inv[row[i]][col[j]] - inv[row[i]][col[k]] * inv[row[k]][col[j]]
          }
        }
        inv[row[i]][col [k]] = -inv[row[i]][col[k]] * inv[row[k]][col[k]]
      }
    }
  }
  // end main reduction loop

  // unscramble rows
  for j in 0..<n {
    for i in 0..<n {
      temp[col[i]] = inv[row[i]][j]
    }
    for i in 0..<n {
      inv[i][j] = temp[i]
    }
  }
  // unscramble columns
  for i in 0..<n {
    for j in 0..<n {
      temp[row[j]] = inv[i][col[j]]
    }
    for j in 0..<n {
      inv[i][j] = temp[j]
    }
  }
  return inv
} // end matinv

// 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
					
print("matrix2.swift running")

var x = [[10.0, 9.0, 8.0], [3.0, 2.0, 1.0]]
print("x= \(x)")
var y = [[1.0, 2.0, 3.0], [4.0, 5.0, 6.0]]
print("y= \(y)")
var z = [[1.0, 2.0, 3.0], [5.0, 6.0, 8.0], [9.0, 11.0,14.0]]
print("z= \(z)")

var sum = x + y
print("sum=x+y= \(sum)")
sum = matadd(x, y)
print("matadd(x,y)= \(sum)")
var inv = matinv(z)
print("matinv(z)= \(inv)")
var det = matdet(z)
print("matdet(z)= \(det)")

print("matrix2.swift finished")
// end matrix2.swift