C Based Reduce Row Echelon Code

Originally Posted on November 7, 2017 by hsauro

I recently needed some code to compute the reduced row echelon of a matrix. Applications such as Matlab, Mathematics, sympy and R support this functionality out of the box. Libraries such as LAPACK do not, including the linear algebra package in Python. The linear algebra package in Python has some surprising omissions which appear to be by design.

Here I include a C based function that will compute the reduced row echelon. It uses partial pivoting to prevent numerical stability but I don’t know how it would fare when confronted with a large matrix (i.e > 500 rows or columns).

It uses a simple matrix type called TMATRIX. It shouldn’t be difficult to modify code if you use a different matrix structure.

To call rref use:

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

The 1e-9 argument is the setting that decides whether a value is zero or not. Numbers below tolerance are considered zeros.

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155

#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); }

The output matrix givne the input should be:

1 2 3 4 5 6

1.000000  0.000000  -2.000000  3.000000  0.000000  -32.666667 0.000000  1.000000  -2.000000  2.000000  0.000000  -9.000000 0.000000  0.000000  0.000000  0.000000  1.000000  5.500000