#include <stdio.h> /* for input/output */
#define NumberOfEquations 4

main() {

  /* VARIABLE DECLARATION */

  /* Data */
  float determinant, sum;
  float matrices[NumberOfEquations][NumberOfEquations+1] = {
    /* A */                                   /* B */
    {  8,  -4,   0,  -2,        6},
    { -4,  12,  -3,   0,        0},
    {  0,  -3,   4,  -1,      -12},
    { -2,   0,  -1,   7,        0}
  };

  /* Indices */
  int substitutionIndex, columnIndex, rowIndex;


  /* INTRODUCTION */

  printf("Lab Session 8, Question 1\n"); /* Display program header */

  /* GAUSS ELIMINATION */

  determinant = 1;
  for (substitutionIndex = 0; substitutionIndex < NumberOfEquations; substitutionIndex++) {
    /* *+ Find maximum pivot in column i then permute rows i and the
       row of the maxiumum pivot */
    determinant *= matrices[substitutionIndex][substitutionIndex];
    for (columnIndex = substitutionIndex+1; columnIndex <= NumberOfEquations; columnIndex++) {
      matrices[substitutionIndex][columnIndex] /= matrices[substitutionIndex][substitutionIndex];
      for (rowIndex = substitutionIndex+1; rowIndex < NumberOfEquations; rowIndex++) {
        matrices[rowIndex][columnIndex] -= 
          matrices[rowIndex][substitutionIndex] * matrices[substitutionIndex][columnIndex];
      }
    }
  }


  /* BACK-SUBSTITUTION */

  for (rowIndex = NumberOfEquations - 2; rowIndex >= 0; rowIndex--) {
    sum = 0;
    for (columnIndex = rowIndex+1; columnIndex < NumberOfEquations; columnIndex++) {
      sum += matrices[rowIndex][columnIndex] * matrices[columnIndex][NumberOfEquations];
    }
    matrices[rowIndex][NumberOfEquations] -= sum;
  }

  for (rowIndex = 0; rowIndex < NumberOfEquations; rowIndex++) {
    printf ("%f\n", matrices[rowIndex][NumberOfEquations]);
  }

}