// test_matinv.swift
//                   matinv.swift should be imported, included below

#if os(OSX) || os(iOS) // for libraries, portable for OSX and Linux
import Foundation  // or Playground
#elseif os(Linux)
import Glibc
#endif 

print("test_matinv.swift running")
var maxerr = 0.0
let n = 5
let m1 = [[1.0, 2.0, 3.0, 4.0, 5.0], 
          [2.0, 3.0, 4.0, 5.0, 1.0],
	  [3.0, 4.0, 5.0, 1.0, 2.0],
	  [4.0, 5.0, 1.0, 2.0, 3.0],
	  [5.0, 4.0, 3.0, 2.0, 1.0]]
var invm1 = [[Double]](repeating:[Double](repeating:0.0,count:n),count:n)
var invinvm1 = [[Double]](repeating:[Double](repeating:0.0,count:n),count:n)

print("m1=\(m1) \n")
invm1 = matinv(m1)
print("matinv(m1)=\(invm1) \n")

invinvm1 = matinv(invm1)
print("matinv(matinv(m1)) close m1=\(invinvm1) \n")

print("check matinv(matinv(m1))-m1 for error")
for i in 0..<n {
  for j in 0..<n {
    maxerr = max(maxerr, abs(m1[i][j]-invinvm1[i][j]))
  }
}
print("maxerr=\(maxerr) \n")
print("test_matinv.swift finished")

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

// end test_matinv.swift