open Impact

let c = [| [| 2.; 1.; -1.; 8.|];
           [| -3.; -1.; 2.; -11.|]; 
           [| -2.; 1.; 2.; -3.|] |] 

let d = [| [| 4.; 8.; 4.|]; 
           [| -3.; -1.; 7.|] |]

let printMatrix(matrix : float[][]) =
    for rowNum in 0 .. matrix.Length-1 do
        let row = matrix.[rowNum]
        for colNum in 0 .. row.Length-1 do
            printf "  %10f" row.[colNum]
        printfn ""
    printfn ""
 
let scalarMult (vector : float [], multiplier : float) =
    [| for i in 0 .. vector.Length-1 -> vector.[i] * multiplier |]
 
let vectorAdd(v1 : float [], v2 : float []) =
    [| for i in 0 .. v1.Length-1 -> v1.[i] + v2.[i]|]

printMatrix c
let scaleRow(m : float [][], d : int) =
    m.[d] <- scalarMult(m.[d], 1.0/m.[d].[d])

let eliminateRow(m : float [][], r: int, d : int) =
    m.[r] <- vectorAdd(scalarMult(c.[d], -c.[r].[d]), m.[r])

let solve3matrix(m : float [][]) = 
 scaleRow(m,0)
 eliminateRow(m,1,0)
 eliminateRow(m,2,0)

 scaleRow(m,1)
 eliminateRow(m,0,1)
 eliminateRow(m,2,1)

 scaleRow(m,2)
 eliminateRow(m,0,2)
 eliminateRow(m,1,2)

let clearCol(m : float [][], d : int) =
  scaleRow(m,d)
  for rowNum in 0 .. m.Length-1 do
    if not (rowNum = d) then
      eliminateRow(m,rowNum,d)

let solve(m : float [][]) =
  for diagonal in 0 .. m.Length-1 do
    clearCol(m,diagonal)  

solve(c)
printMatrix c