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