/* Dynamic 2D Number Array with Gauss Elimination Support; version 1.0 */
/* (c) copyright Ian Hickson 2000, distributed under the GNU GPL */

#include <math.h> /* fabs */
#include "dynamic-2d-arrays.h" /* array2d stuff */

/* Implements a float dynamic square array with gauss elimination */

typedef float geNumber;

array2d* geAllocArray(int equations) {
  array2d* array;
  array = allocArray2D(equations+1, equations, sizeof(geNumber));
  return array;
}

void geSet(array2d* array, int x, int y, geNumber value) {
  setArray2D(array, x, y, &value);
}

geNumber geGet(array2d* array, int x, int y) {
  geNumber value;
  value = *(geNumber*)getArray2D(array, x, y);
  return value;
}

void gePrint(array2d* array, int dp) {
  int x, y;
  for(y = 0; y < array->height; y++) {
    for(x = 0; x < array->width; x++) {
      printf("%6.*f", dp, geGet(array, x, y));
    }
    printf("\n");
  }
}

int geGetMaxPivot(array2d* array, int pivot) {
  int selectedRow = pivot;
  int y;
  geNumber maximumPivotValue = 0;
  geNumber value;
  for (y = pivot; y < array->height; y++) {
    value = fabs(*(geNumber*)getArray2D(array, pivot, y));
    if (value > maximumPivotValue) {
      maximumPivotValue = value;
      selectedRow = pivot;
    }
  }
  return selectedRow;
}

array2d* geGaussElimination(array2d* originalArray) {
  int substitutionIndex, x, y;
  geNumber determinant = 1;
  geNumber value, sum;
  array2d* answers;
  array2d* array;
  array = cloneArray2D(originalArray);
  answers = allocArray2D(1, array->height, array->elementSize);
  for (substitutionIndex = 0; substitutionIndex < array->height; substitutionIndex++) {
    /* Find maximum pivot in column i then premute row i and the row of the maximum pivot: */
    y = geGetMaxPivot(array, substitutionIndex);
    permuteRowsArray2D(array, y, substitutionIndex);
    determinant *= pow(-1, y-substitutionIndex) * geGet(array, substitutionIndex, substitutionIndex);
    for (x = substitutionIndex+1; x < array->width; x++) {
      geSet(array, x, substitutionIndex, 
            geGet(array, x, substitutionIndex) / geGet(array, substitutionIndex, substitutionIndex)); 
      for (y = substitutionIndex+1; y < array->height; y++) {
        geSet(array, x, y, 
              geGet(array, x, y) - (geGet(array, substitutionIndex, y) * geGet(array, x, substitutionIndex))); 
      }
    }
  }
  /* Backsubstitution */
  for (substitutionIndex = array->height-2; substitutionIndex >= 0; substitutionIndex--) {
    sum = 0;
    for (x = substitutionIndex+1; x < array->width-1; x++) {
      sum += geGet(array, x, substitutionIndex) * geGet(array, array->width-1, x);
    }
    geSet(array, array->width-1, substitutionIndex, (geGet(array, array->width-1, substitutionIndex) - sum));
  }
  /* Answer Collection */
  for (y = 0; y < array->height; y++) {
    geSet(answers, 0, y, geGet(array, array->width-1, y));
  }
  freeArray2D(array);
  return answers;
}

/* end */