/*
{m:nat, n:nat}
record <'a> matrix(m,n) {
  row: int(m);
  col: int(n);
  data[m][n]: 'a
}
*/

float fabs (f: float) {
  if (f >. 0.0) return f; else return -.f;
}

('a){m:nat,n:nat,i:nat,j:nat | i < m, j < m}
unit
rowSwap(data[m][n]: 'a, i:int(i), j:int(j)) {
  var: temp[]: 'a;;
  temp = data[i];
  data[i] = data[j];
  data[j] = temp;
}

{n:nat,i:nat | i < n}
unit
norm(r[n]: float, n:int(n), i:int(i)) {
  var: float x;;

  x = r[i]; r[i] = 1.0; i = i + 1;

  invariant: [i:nat] (i: int(i))
  while (i < n) { r[i] = r[i] /. x; i = i + 1;}
}

{n:nat, i:nat | i < n}
unit
rowElim(r[n]:float, s[n]:float, n:int(n), i:int(i)) {
  var: float x;;

  x = s[i]; s[i] = 0.0; i = i + 1;

  invariant: [i:nat] (i: int(i))
  while (i < n) { s[i] = s[i] -. x *. r[i]; i = i + 1;}
}

{m:nat, n:nat, i:nat | m > i, n > i}
int[0, m)
rowMax (data[m][n]: float, m: int(m), i: int(i)) {
  var: nat j; float x, y;
       max: int[0, m);;

  max = i; j = i + 1;
  x = fabs(data[i][i]); /* fabs: absolute value function for float */
  
  invariant: [j:nat] (j: int(j))
  while (j < m) {
    y = fabs(data[j][i]);
    if (y >. x) { x = y; max = j; }
    j = j + 1;
  }
  return max;
}

{n:nat | n > 0}
unit gauss (mat: <float> matrix(n,n+1)) {
  var: float data[][n+1]; nat i, j, max, n;;

  data = mat.data; n = mat.row;
  for (i = 0; i < n; i = i + 1) {
    max = rowMax(data, n, i);
    norm (data[max], n+1, i);
    rowSwap(data, i, max);
    for (j = i + 1; j < n; j = j + 1) {
      rowElim(data[i], data[j], n+1, i);
    }
  }
}

{n:nat}
unit print_array (size: int(n), a[n]: float) {
  var: nat i;;
  for (i = 0; i < size; i = i + 1) {
      print_float (a[i]); print_string ("\t");
    }
    print_newline ();
}

unit print_matrix (a: <float> matrix) {
  var: nat i;;
  for (i = 0; i < a.row; i = i + 1) {
      print_array (a.col, a.data[i]);
    }
}

unit main () {
  var: a: <float> matrix (3, 4);;

  a = matrix { row = 3; col = 4; data = alloc (3, alloc (4, 0.0)) };

  a.data[0][0] = 1.0; a.data[0][1] = 2.0;
  a.data[0][2] = -.3.0; a.data[0][3] = -.4.0;

  a.data[1] = alloc (4, 0.0);
  a.data[1][0] = -.2.0; a.data[1][1] = 3.0;
  a.data[1][2] = 1.0; a.data[1][3] = 7.0;

  a.data[2] = alloc (4, 0.0);
  a.data[2][0] = 3.0; a.data[2][1] = -.1.0;
  a.data[2][2] = 2.0; a.data[2][3] = 7.0;

  print_matrix (a);
  print_newline ();
  gauss (a);
  print_matrix (a);
}