Jump to content
  • Advertisement
Sign in to follow this  
EvanDavis

Gauss-Jordan Problems

This topic is 2440 days old which is more than the 365 day threshold we allow for new replies. Please post a new topic.

If you intended to correct an error in the post then please contact us.

Recommended Posts

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

Share this post


Link to post
Share on other sites
Advertisement
Sign in to follow this  

  • Advertisement
×

Important Information

By using GameDev.net, you agree to our community Guidelines, Terms of Use, and Privacy Policy.

We are the game development community.

Whether you are an indie, hobbyist, AAA developer, or just trying to learn, GameDev.net is the place for you to learn, share, and connect with the games industry. Learn more About Us or sign up!

Sign me up!