#include<stdio.h>

#include<stdlib.h>

#include<math.h>

typedef struct TMatrix {

double **m;

int nRows; int nCols;

} TMATRIX;

TMATRIX createMatrix (int nRows, int nCols) {

TMATRIX oMatrix;

oMatrix.nRows = nRows;

oMatrix.nCols = nCols;

oMatrix.m = (double **) malloc (sizeof (double*) * sizeof (double*));

for (int i = 0; i < nRows; i++)

oMatrix.m[i] = (double*) malloc (sizeof (double)*nCols);

return oMatrix;

}

// Remove values less than tolerance

void cleanUpMatrix (TMATRIX oMatrix, double tolerance) {

// Removes all numbers close to zero, i.e between -tol and +tol

for (int i = 0; i < oMatrix.nRows; i++)

for (int j = 0; j < oMatrix.nCols; j++)

if (fabs (oMatrix.m[i][j]) < tolerance)

oMatrix.m[i][j] = 0;

}

void exchangeRows (TMATRIX oMatrix, int r1, int r2) {

double t = 0;

for (int i = 0; i < oMatrix.nCols; i++) {

t = oMatrix.m[r1][i];

oMatrix.m[r1][i] = oMatrix.m[r2][i];

oMatrix.m[r2][i] = t;

}

}

// Mainly fo rdeugging, helps spott out of range errors when indexing elements

double getValue (TMATRIX oMatrix, int i, int j) {

if ((i < 0) || (j < 0)) {

printf ("Error in indexing\n");

getchar ();

exit (0);

}

if ((i >= oMatrix.nRows) || (j >= oMatrix.nCols)) {

printf ("Error in indexing: %d, %d\n", i, j);

getchar ();

exit (0);

}

return oMatrix.m[i][j];

}

void printMatrix (TMATRIX oMatrix) {

for (int i = 0; i < oMatrix.nRows; i++) {

for (int j = 0; j < oMatrix.nCols; j++)

printf ("%f ", oMatrix.m[i][j]);

printf ("\n");

}

}

TMATRIX rref (TMATRIX oMatrix, double tolerance)

{

int currentRow; double factor;

TMATRIX oEchelon = createMatrix (oMatrix.nRows, oMatrix.nCols);

// Make a copy and work on that.

for (int i = 0; i < oMatrix.nRows; i++)

for (int j = 0; j < oMatrix.nCols; j++)

oEchelon.m[i][j] = oMatrix.m[i][j];

int Arow = 0; int Acol = 0;

while ((Arow < oEchelon.nRows) && (Acol < oEchelon.nCols)) {

// locate a nonzero column

if (fabs (getValue (oEchelon, Arow, Acol) < tolerance)) {

// If the entry is zero work our way down the matrix

// looking for a nonzero entry, when found, swap it for Arow

currentRow = Arow;

do {

// next row

currentRow++;

// Have we reached the end of the rows but we've still got columns left to scan?

if ((currentRow >= oEchelon.nRows) && (Acol <= oEchelon.nCols)) {

// reset row counter back to where it was and try next column

currentRow = Arow; Acol++;

}

// If we've scanned the whole matrix, then lets get out...

if (currentRow >= oEchelon.nRows) {

cleanUpMatrix (oEchelon, tolerance);

return oEchelon;

}

} while (fabs (getValue (oEchelon, currentRow, Acol)) < tolerance);

// We've found a nonzero row entry so swap it with 'Arow' which did have a zero as its entry

exchangeRows (oEchelon, Arow, currentRow);

}

// Arow now holds the row of interest }

factor = 1.0 / getValue (oEchelon, Arow, Acol);

// reduce all the entries along the column by the factor

for (int i = Acol; i < oEchelon.nCols; i++)

oEchelon.m[Arow][i] = getValue (oEchelon, Arow, i) * factor;

// now eliminate all entries above and below Arow, this generates the reduced form

for (int i = 0; i < oEchelon.nRows; i++) {

// miss out Arow itself

if ((i != Arow) && (fabs (getValue (oEchelon, i, Acol)) > tolerance)) {

factor = getValue (oEchelon, i, Acol);

// work your way along the column doing the same operation

for (int j = Acol; j < oEchelon.nCols; j++) {

oEchelon.m[i][j] = getValue (oEchelon, i, j) - factor * getValue (oEchelon, Arow, j);

}

}

}

Arow++; Acol++;

}

cleanUpMatrix (oEchelon, tolerance);

return oEchelon;

}

int main (int argc, const char* argv[]) {

TMATRIX oEchelon;

printf ("\nTesting rref\n\n");

TMatrix oMatrix = createMatrix (3, 6);

oMatrix.m[0][0] = 0; oMatrix.m[0][1] = 3; oMatrix.m[0][2] = -6; oMatrix.m[0][3] = 6; oMatrix.m[0][4] = 4; oMatrix.m[0][5] = -5;

oMatrix.m[1][0] = 3; oMatrix.m[1][1] = -7; oMatrix.m[1][2] = 8; oMatrix.m[1][3] = -5; oMatrix.m[1][4] = 8; oMatrix.m[1][5] = 9;

oMatrix.m[2][0] = 3; oMatrix.m[2][1] = -9; oMatrix.m[2][2] = 12; oMatrix.m[2][3] = -9; oMatrix.m[2][4] = 6; oMatrix.m[2][5] = 16;

oEchelon = rref (oMatrix, 1e-6);

for (int j = 0; j < 20; j++) {

printMatrix (oEchelon);

printf ("\n");

}

getchar ();

exit (0);

}