Gaussian elimination

[lekfanda]

Hi guys,

Ok so the excercise I'm having trouble with is a Gaussian Elimination.
So I'm given a matrix and the vector in respectively 2 txt files: A.txt and b.txt for the equation A * x = b .
So I scanned and stored the matrix and vector into A[500][500] and b[500]. I assure you everything checks up until the Gaussian Elimination loop below. I don't know what's wrong but when I check the changed matrix I can see the diagonal but most of the numbers below the Diagonal are only super close to 0 and not exactly 0.

 for (int i = 0; i < p-1; i++)
	{
		for (int j = i+1; j < p; j++)
		{
			double c = A[j][i] / A[i][i];
			for (int k = 0;k < p;k++) {
				
				A[j][k] = A[j][k] - c * A[i][k];
			}
		}
	}

ofstream pivot5("pivot5.txt");
	for (int i = 0;	i < p;	i++) {
		for (int j = 0;  j < p;  j++)
		{
			pivot5 << A[i][j] << " ";
		}
		pivot5 << "\n";
	}

Im posting the output below.
My assignment is due in 20 hours so I would appreciate any quick tip.

Link to the output file "pivot5.txt": https://jpst.it/17K7X (if you copy paste it into a larger txt editor its easier to see the diagonal and the pattern)

[jonnin]

it is normal to get near but not exactly zero. This is normal roundoff and your code is probably fine if it gives approximately the correct answer.

There are numerical methods to improve the answer and get closer if you need this for a real application. If this is classroom work, just let it be. I would not expect a brute force GE solver to provide the exact answer for anything this large -- it can do it on smaller matrices like 20x20 or smaller, but 500x500 is asking a lot from this approach.

[lastchance]

(1) Test your program with a small matrix first - say 3*3. It doesn't make sense to start with 500*500.

(2) Try your program on the matrix

0 1 0
3 4 5
6 7 8

There will be a problem, even though this matrix is non-singular. Hint: think what you are doing with the diagonal element.

(3) If you are genuinely doing gaussian elimination then you should simultaneously be doing the same pivoting on the vector b; I can see no evidence of it, nor of the back-substitution afterwards to find x.

(4) Your ultimate test is whether you can premultiply your solution x by A and get b.

(5) With rows of length 500, round-off errors will accumulate rapidly. You could add a line after line 9 deliberately setting A[j][i]=0. Also, since you will already have cleared all columns before the ith, you could start the loop on line 6 at k=i, rather than k=0. Overall, this will half the amount of computation, as well as reducing round-off error.

[lekfanda]

Ok, so I checked the assignment once again and it was actually written that for this excercise the diagonal values may never reach 0 so that part is fine.

(3) If you are genuinely doing gaussian elimination then you should simultaneously be doing the same pivoting on the vector b; I can see no evidence of it, nor of the back-substitution afterwards to find x.

I also just did that so now I have done the pivoting to the Matrix A and the results Vector b (Output for the pivoted b: https://jpst.it/17LQc ) through this code:

	for (int i = 0; i < p-1; i++)
	{
		for (int j = i+1; j < p; j++)
		{
			double c = A[j][i] / A[i][i];
			for (int k = 0;k < p;k++) {
				
				A[j][k] = A[j][k] - c * A[i][k];
				
			}
			b[j] = b[j] - c*b[i];
		}
	}

and it checks out (I even manually checked the first 2 lines), but I'm having problem with the solution "x" for which I used this code:

	double x[p];

	for (int i = p-1; i > 0; i--)
	{
		x[i] = A[i][p];
		for (int j = i+1; j < p; j++)
		{
			if (j != i) x[i] = b[i] - A[i][j] * x[j];
			x[i] = x[i] / A[i][i];
		}
	}

When I multiply x by A to test it I don't get "b" but some really big values instead.

ofstream xtxt("x.txt");
	for (int i = 0; i < p; i++)
	{
		xtxt << x[i] <<endl;
	}
	double c[p];
	for (int i = 0; i < k; i++) {
		c[i] = 0;
		for (int j = 0; j < l; j++)
		{
			c[i] += A[i][j] * x[j];
		}
	}
	ofstream ctxt("c.txt");
	for (int i = 0; i < p; i++)
	{
		ctxt << c[i] << endl;
	}

x.txt: https://jpst.it/17LRX
c.txt: https://jpst.it/17LS3
b.txt:https://jpst.it/17LS7

[lastchance]

I hope that you have been doing Gaussian elimination on a COPY of A and b, because you change the original matrices.

For the back-substitution:

   // Back-substitute
   for ( int i = p - 1; i >= 0; i-- )
   {
      x[i] = b[i];
      for ( int k = i + 1; k < p; k++ ) x[i] -= A[i][k] * x[k];
      x[i] /= A[i][i];
   }

In your back-substitution code (the middle sample of your post) on line 5, A[i][p] would be going outside array bounds, line 8 should be cumulating changes to x[i], rather than replacing them, and what is currently line 9 should be after the j loop. You should think how you do this by hand with a small example.

[lekfanda]

I hope that you have been doing Gaussian elimination on a COPY of A and b, because you change the original matrices.

For the back-substitution:
1
2
3
4
5
6
7
// Back-substitute
for ( int i = p - 1; i >= 0; i-- )
{
x[i] = b[i];
for ( int k = i + 1; k < p; k++ ) x[i] -= A[i][k] * x[k];
x[i] /= A[i][i];
}


In your back-substitution code (the middle sample of your post) on line 5, A[i][p] would be going outside array bounds, line 8 should be cumulating changes to x[i], rather than replacing them, and what is currently line 9 should be after the j loop. You should think how you do this by hand with a small example.

Yep, you were right about the copy of A and b: that was the reason why the test multiplication wouldn't give me the 'b' vector. I also fixed the back-substitution according to your instructions and everything worked fine. Thanks for another solved one ;)