• Create Account

FREE SOFTWARE GIVEAWAY

We have 4 x Pro Licences (valued at \$59 each) for 2d modular animation software Spriter to give away in this Thursday's GDNet Direct email newsletter.

Read more in this forum topic or make sure you're signed up (from the right-hand sidebar on the homepage) and read Thursday's newsletter to get in the running!

# Gauss-Jordan Problems

Old topic!
Guest, the last post of this topic is over 60 days old and at this point you may not reply in this topic. If you wish to continue this conversation start a new topic.

### #1cptroot  Members   -  Reputation: 100

Like
0Likes
Like

Posted 13 December 2011 - 12:25 AM

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


### #2inprazia  Members   -  Reputation: 107

Like
0Likes
Like

Posted 16 December 2011 - 08:17 AM

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

Old topic!
Guest, the last post of this topic is over 60 days old and at this point you may not reply in this topic. If you wish to continue this conversation start a new topic.

PARTNERS