I've been working on some code for finding the inverse of a matrix, but I keep running into problems with my implementation of Gauss-Jordan Elimination.
Ordinarily, Gauss-Jordan Elimination finds the inverse of a matrix using row operations on a matrix, but I keep ending up with blank rows. Whenever I try to find the inverse of a matrix that has more than 4 rows, I get empty rows that make my algorithm think the matrix is singular, even when Wolfram-Alpha tells me it's not.

Here's the code, simplified. Can anyone help find the problem.



Matrix inverse() {
if (length[0] != length[1]) throw new Exception("Non-Square Matrix");
Matrix augmented = this.augment(identity(length[1]));
double eps = .00000000001;

int h = length[0];
int w = length[0] * 2;
double c;
foreach (y; 0..h) {
int maxrow = y;
foreach (i; y+1..h) { // Find max pivot
if (to!double(abs(augmented[i, y])) > to!double(abs(augmented[maxrow, y]))) {
maxrow = i;
}
}
int maxrow = y + 1;
double[] temp;
foreach (j; 0..w) {
temp ~= augmented[y, j];
}
foreach (j; 0..w) {
augmented[y, j] = augmented[maxrow, j];
}
foreach (j; 0..w) {
augmented[maxrow, j] = temp[j];
}

if (to!double(abs(augmented[y, y])) <= eps) // Singular?
return Matrix.init;
foreach (i; y+1..h) { // Eliminate column y
c = augmented[i, y] / augmented[y, y];
foreach (x; y..w) {
augmented[i, x] -= augmented[y, x] * c;
if (abs(augmented[i, x]) <= eps) augmented[i, x] = 0;
}
}
}
for (int y = h-1; y > -1; y--) { // Backsubstitute
c = augmented[y, y];
foreach (i; 0..y) {
for (int x = w - 1; x > y - 1; x--) {
augmented[i, x] -= augmented[y, x] * augmented[i, y] / c;
if (abs(augmented[i, x]) <= eps) augmented[i, x] = 0;
}
}
augmented[y, y] /= c;
foreach (x; h..w) { // Normalize row y
augmented[y, x] /= c;
}
}
Matrix result = Matrix(h, h);
foreach (i; 0..h) {
foreach (j; 0..h) {
result[i, j] = augmented[i, j + h];
}
}
return result;
}


At first glance it seems that you should eliminate that peace of code: int maxrow = y + 1;