// matrix.swift           should be importable, test included at bottom
#if os(OSX) || os(iOS)
import Darwin
#elseif os(Linux)
import Glibc
#endif

print("matrix.swift running, should be an import, test included")

// vector functions
// absolute value of a vector
func vecabs(_ numbers: [Double]) -> Double {
  var total: Double = 0
  for number in numbers {
    total += abs(number)
  }
  return total
} // end vecabs

// length of a vector, sum of squares
func veclen(_ numbers: [Double]) -> Double {
  var total: Double = 0
  for number in numbers {
    total += number*number
  }
  return total
} // end veclen

// vector difference errors average, RMS, max
func vecerr(_ a: [Double], _ b: [Double]) {
  let n = a.count
  var sum = 0.0
  var sumsq = 0.0
  var adiff = 0.0
  var max = abs(a[0]-b[0])
  for i in 0..<n {
    adiff = abs(a[i]-b[i])
    sum = sum + adiff
    sumsq = sumsq + adiff*adiff
    if adiff>max {
      max = adiff
    }
  }
  sumsq = sqrt(sumsq/Double(n))
  adiff = sum/Double(n)
  print("errors average=\(adiff), RMS=\(sumsq), max=\(max)")
  return
} // end vecerr

// vector dot product 
func vecdot(_ a: [Double], _ b: [Double]) -> Double {
  let n = a.count
  var dot = 0.0
  for i in 0..<n {
    dot = dot + a[i]*b[i] 
  }
  return dot
} // end vecdot

// vector cross product 
func veccross(_ x: [Double], _ y: [Double]) -> [Double] {
  let n = x.count
  var prod = [Double](repeating:0.0,count:n)
  var k = 0
  var sign = 1.0
  for m in 0..<n {
    for j in 0..<n {
      k = j + m
      if k>=n { k = k-n }
      prod[m] = prod[m] + x[j]+y[k]
    }
    prod[m] = sign*prod[m]
    sign = -sign
  }
  return prod
} // end veccross

// vector sum 
func vecadd(_ a: [Double], _ b: [Double]) -> [Double] {
  let n = a.count
  var sum = [Double](repeating:0.0,count:n)
  for i in 0..<n {
    sum[i] = a[i]+b[i] 
  }
  return sum
} // end vecadd

// vector difference 
func vecsub(_ a: [Double], _ b: [Double]) -> [Double] {
  let n = a.count
  var dif = [Double](repeating:0.0,count:n)
  for i in 0..<n {
    dif[i] = a[i]-b[i] 
  }
  return dif
} // end vecsub

// matrix functions
// make a matrix   n by n  zero values
func matmake(_ n: Int) -> [[Double]] {
  let mat = [[Double]](repeating:[Double](repeating:0.0,count:n),count:n)
  return mat
} // end matmk

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

// add two n by m matrices
func mataddv(_ a: [[Double]], _ b: [[Double]]) -> [[Double]] {
  let nia = a.count
  // print("mataddv nia = \(nia)")
  let nja = a[0].count
  // print("mataddv nja = \(nja)")
  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]
    }
  }
  // print("mataddv  nia = \(nia)  nja = \(nja) \n")
  return sum
} // end mataddv

// subtract matrix
func matsub(_ n: Int,_ a: [[Double]], _ b: [[Double]]) -> [[Double]] {
  var dif = [[Double]](repeating:[Double](repeating:0.0,count:n),count:n)
  for i in 0..<n {
    for j in 0..<n {
      dif[i][j] = a[i][j]-b[i][j]
    }
  }
  return dif
} // end matsub

// multiply matrix times matrix, may not be square
func matmulv(_ a: [[Double]], _ b: [[Double]]) -> [[Double]] {
  let ni = a.count
  let nk = a[0].count
  let nm = b.count
  let nj = b[0].count
  if nk != nm || ni != nj {
    print("inconsistant matrix sizes a=\(ni),\(nk) b=\(nm),\(nj)")
    return  [[Double]](repeating:[Double](repeating:0.0,count:nk),count:ni)
  }
  var prod = [[Double]](repeating:[Double](repeating:0.0,count:ni),count:nj)
  for i in 0..<ni {
    for j in 0..<nk {
      prod[i][j] = 0.0
      for k in 0..<nk {
        prod[i][j] = prod[i][j] + a[i][k]*b[k][j]
      }
    }
  }
  return prod
} // end matmulv

// 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


// a function to print a matrix
func matprt(_ n: Int, _ mat: [[Double]]) {
  for i in 0..<n {
    for j in 0..<n {
      print("mat[\(i)][\(j)] = \(mat[i][j])")
    }
  }
} // end matprt

// a function to print a non square matrix
func matprtv(_ a: [[Double]]) {
  let nia = a.count
  let nja = a[0].count
  for i in 0..<nia {
    for j in 0..<nja {
      print("mat[\(i)][\(j)] = \(a[i][j])")
    }
  }
} // end matprtv

print("\n test vector functions \n")
var a = [1.0, 2.0, 3.0, 4.0]
var b = [2.0, 4.0, 6.0, 8.0]
print("a = \(a) \n")
print("b = \(b) \n")
var c = vecadd(a, b)
print("c=a+b= \(c) \n")
var d = vecsub(a, b)
print("d=a-b= \(d) \n")
print("veclen(a) = \(veclen(a))")
print("veclen(b) = \(veclen(b))")
print("veclen(c) = \(veclen(c))")
print("veclen(d) = \(veclen(d))")
print("vecabs(d) = \(vecabs(d))")
var ar = [1.0, 2.2, 3.5, 3.4]
vecerr(a, ar)
print(" ")
print("vecdot(a,b) = \(vecdot(a,b))")
var e = veccross(a,b)
print("e=veccross(a,b) = \(e))")

print("\n test matrix functions \n")

let A = [[1.0, 2.0, 3.0],[5.0, 4.0, 6.0], [7.0, 9.0, 8.0]]
print("A = \(A)")
matprt(3, A)
print("A[2][2] = \(A[2][2]) \n")
print("A.count = \(A.count)")
print("A[0].count = \(A[0].count)")

var Av = [[1.0, 2.0],[4.0, 5.0], [7.0, 8.0]]
var Bv = [[1.0, 2.0, 3.0], [4.0, 5.0, 6.0]]
print("\n Av.count = \(Av.count)")
print("\n Av[0].count = \(Av[0].count)")
print("Av=")
matprtv(Av)
print("Av[1][1] = \(Av[1][1]) \n")
print("Bv=")
matprtv(Bv)

var n = 3
for i in 0..<n {
  for j in 0..<n {
    print("i=\(i) j=\(j)  A=\(A[i][j])")
  }
}


print("\n now loops to initialize")
var B = [[Double]](repeating:[Double](repeating:0.0,count:3),count:3)
// zeroed 3 by 3 matrix for  Double 
for i in 0..<n {
  for j in 0..<n {
    B[i][j] = Double(i*n+j)
  }
}
for i in 0..<n {
  for j in 0..<n {
    print("i=\(i) j=\(j)  B=\(B[i][j])")
  }
}

print("\nmatadd  add two matrix")
var C = matmake(3)
C = matadd(3, A, B)
print("matprt  print sum of A+B")
matprt(3, C)

print("\nmataddv  add two matrix")
var Cv = mataddv(A, B)
print("\nmatprt  print C = sumv of A+B")
matprtv(Cv)

print("\nmatsub  subtract matrix")
C = matsub(3, A, B)
print("matprt  print C = difference of A-B")
matprt(3, C)

C = matmulv(A, B)
print("\nmatprtv  print C = matmulv product of A*B")
matprt(3, C)

var D = matinv(A)
print("\nmatprtv  print D = matinv inverse of A")
matprtv(D)

var ID = matmulv(A, D)
print("\nmatprtv  print ID = A * D = identity")
matprtv(ID)

print("\n matrix.swift finished")