Jump to content

  • Log In with Google      Sign In   
  • Create Account


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.

  • You cannot reply to this topic
1 reply to this topic

#1 cptroot   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;
}


Sponsor:

#2 inprazia   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