Started by Mar 05 2004 05:38 AM

,
27 replies to this topic

Posted 05 March 2004 - 05:38 AM

At this moment I am coding a 4x4 matrix class. I am trying to write an algorithm for doing the inverse using the Gaussian Elimination method. I have found a few ok resources on the internet about this method but nothing that really helps me 100%.
I was wondering if anyone has coded a 4x4 matrix class and used gaussian elimination to find the inverse and if so if they would be willing to discuss it with me. Below are some of the links I have found on the subject and a pretty good link on Gauss-Jordan elimination but I have one major problem with the Gauss-Jordan elimination.
http://mathworld.wolfram.com/GaussianElimination.html
http://www.mcraefamily.com/mathhelp/MatrixReducedRowEchelonForm.htm
http://www.rism.com/LinAlg/LAsix.htm
http://www.rism.com/LinAlg/LAtwo.htm
Thank you for all your help in advance.
optics
[edited by - optics on March 7, 2004 7:29:13 AM]

Posted 05 March 2004 - 08:55 AM

quote:Original post by optics

... but I have one major problem with the Gauss-Jordan elimination.

If you mention what the problem is we might be able to help you better.

Posted 05 March 2004 - 09:01 AM

We have to solve a 9x9 matrix by next Friday for class, I''ll tell you how it goes.

Posted 05 March 2004 - 09:43 AM

The problem i have with Gauss-Jordan(a simplified version of Gaussian) is that you can end up with a fourth row of 0 0 0 0.

It is hard to mulitply anything by 0 and get 1 for the identity matrix. So this poses a problem for me unless matrices that end up with a 4th row 0 0 0 0 have no determinant.

The problem that I have in coding the algorithm for Gaussian elimination is that everything I find on the internet is that they the pages are mostly theory and not a lot of problem solving. If I could see the Gaussian elimination done on a mock up problem then it would not be that difficult to code. If someone can show me a resource that maybe i missed that works out a sample problem and explains the steps then I could code this algorithm. After I am done then I would be glad to post the source here and also the explanation so that people after me can understand what is going on.

ManBot1: There are NxN matrices classes out there on the internet that would be able to solve your problem. My problem is specific on a 4x4 using Gaussian Elimination. Unless your 9x9 problem will be using the Gaussian Elimination method. If you know a lot about this method or have a resource on it please let me know.

It is hard to mulitply anything by 0 and get 1 for the identity matrix. So this poses a problem for me unless matrices that end up with a 4th row 0 0 0 0 have no determinant.

The problem that I have in coding the algorithm for Gaussian elimination is that everything I find on the internet is that they the pages are mostly theory and not a lot of problem solving. If I could see the Gaussian elimination done on a mock up problem then it would not be that difficult to code. If someone can show me a resource that maybe i missed that works out a sample problem and explains the steps then I could code this algorithm. After I am done then I would be glad to post the source here and also the explanation so that people after me can understand what is going on.

ManBot1: There are NxN matrices classes out there on the internet that would be able to solve your problem. My problem is specific on a 4x4 using Gaussian Elimination. Unless your 9x9 problem will be using the Gaussian Elimination method. If you know a lot about this method or have a resource on it please let me know.

Posted 05 March 2004 - 09:56 AM

If you end with a row of all 0s then the matrix is non-invertible. So that really isn''t an issue.

Posted 06 March 2004 - 01:12 AM

Are you trying to invert an arbitary matrix?

Note that transformation matrices (such as rotation matrices)

are orthogonal and thus can be inverted by simply transposing

them.

In most 3D applications this is what you want so you don''t have

to implement general matrix inversion.

Note that transformation matrices (such as rotation matrices)

are orthogonal and thus can be inverted by simply transposing

them.

In most 3D applications this is what you want so you don''t have

to implement general matrix inversion.

Posted 06 March 2004 - 01:53 AM

Given you''re using a fixed size of 4 x 4 wouldn''t it be just as or similarly efficent to calculate the inverse in the ''classic'' way ? ie:

I mean 4x4 is pretty small, you''d avoid any branching code (mis-predictions etc). You''d be looking at 56 mults + 23 sub/add instructions to complete the inversion, and I''m pretty sure if you played arround a bit on paper you could reduce the number of mults...

(1/det(A)) * TRANSFORM( A )

I mean 4x4 is pretty small, you''d avoid any branching code (mis-predictions etc). You''d be looking at 56 mults + 23 sub/add instructions to complete the inversion, and I''m pretty sure if you played arround a bit on paper you could reduce the number of mults...

Posted 06 March 2004 - 06:44 PM

darookie: The matrix is not of arbitrary size.

I was not aware of that for transformation matrices. That may prove to be the quicker route to go but I also want to do this just in case a matrix inversion is necessary on down the road.

Anonymous: I have heard that Gaussian Elimination is a more efficient way of calculating an inverse matrix, although I have not actually been able to test this theory I can not give an exact answer. We will find out when I finish coding this thing.

i will be using an evaluation tool by intel to discover the truth. http://www.gamasutra.com/features/20000131/barad_pfv.htm

at the end of this document is the tool.

On another note if anyone knows the answer to this particular problem on speed, Gaussian Elimination vs. 1/det(A)* Transform (A) Please share with us. It would probably save me some headaches and time.

Thanks.

[edited by - optics on March 7, 2004 7:32:39 AM]

I was not aware of that for transformation matrices. That may prove to be the quicker route to go but I also want to do this just in case a matrix inversion is necessary on down the road.

Anonymous: I have heard that Gaussian Elimination is a more efficient way of calculating an inverse matrix, although I have not actually been able to test this theory I can not give an exact answer. We will find out when I finish coding this thing.

i will be using an evaluation tool by intel to discover the truth. http://www.gamasutra.com/features/20000131/barad_pfv.htm

at the end of this document is the tool.

On another note if anyone knows the answer to this particular problem on speed, Gaussian Elimination vs. 1/det(A)* Transform (A) Please share with us. It would probably save me some headaches and time.

Thanks.

[edited by - optics on March 7, 2004 7:32:39 AM]

Posted 06 March 2004 - 06:54 PM

Gaussian elimination : less calculations, more memory accesses

(1/det)*adjoint : more calculations, less memory accesses

I can't remember where the turnover point is, but for large matrices the number of calculations required by the (1/det)*adjoint method end up taking more time than the memory accesses for gaussian elimination.

Large matrices : Gaussian elimination is faster

Small matrices : (1/det)*adjoint is faster

I think I remember the (1/det)*adjoint method being faster for 4x4 matrices, but don't quote me on that.

[edited by - joanusdmentia on March 7, 2004 1:55:36 AM]

(1/det)*adjoint : more calculations, less memory accesses

I can't remember where the turnover point is, but for large matrices the number of calculations required by the (1/det)*adjoint method end up taking more time than the memory accesses for gaussian elimination.

Large matrices : Gaussian elimination is faster

Small matrices : (1/det)*adjoint is faster

I think I remember the (1/det)*adjoint method being faster for 4x4 matrices, but don't quote me on that.

[edited by - joanusdmentia on March 7, 2004 1:55:36 AM]

Posted 07 March 2004 - 12:01 AM

Ok folks. Got the hang of the calculations that are needed to perform the Gaussian elimination however I know that I am not getting the correct form. I feel this way because of this site

http://mathworld.wolfram.com/GaussianElimination.html

What is bugging me is the b_{1} ... b_{k} part of the equation.

Suppose we have a 4x4 matrix

| e11 e12 e13 e14 |

| e21 e22 e23 e24 |

| e31 e32 e33 e34 |

| e41 e42 e43 e44 |

and we wanted to use Gaussian Elimination to get the matrix in Triangular form

| e11 e12 e13 e14 |

| 0 ~e22 ~e23 ~e24 |

| 0 0 ~e33 ~e34 |

| 0 0 0 ~e44 |

then solve for x_{1} to x_{k}. I am having a hard time trying to figure out what b_{1}, b_{2}, b_{3}, b_{4} is suppose to be set to. Do we pull these values out of thin air or is there something I overlooked?

Getting closer to the finish line.

[edited by - optics on March 7, 2004 7:34:06 AM]

[edited by - optics on March 7, 2004 10:55:56 AM]

http://mathworld.wolfram.com/GaussianElimination.html

What is bugging me is the b

Suppose we have a 4x4 matrix

| e11 e12 e13 e14 |

| e21 e22 e23 e24 |

| e31 e32 e33 e34 |

| e41 e42 e43 e44 |

and we wanted to use Gaussian Elimination to get the matrix in Triangular form

| e11 e12 e13 e14 |

| 0 ~e22 ~e23 ~e24 |

| 0 0 ~e33 ~e34 |

| 0 0 0 ~e44 |

then solve for x

Getting closer to the finish line.

[edited by - optics on March 7, 2004 7:34:06 AM]

[edited by - optics on March 7, 2004 10:55:56 AM]

Posted 07 March 2004 - 12:02 AM

quote:Original post by optics

darookie: The matrix is of abitrary size( i.e. 4x4 ) but the values are not.

Not wishing to critisise/be unfair etc., but I''m assuming english isn''t your first language - you seem to be mis-understanding ''arbitrary'' - which means something that could be anything (it hasn''t been specified) of a particular class (for example "an arbitrary distance" would be something that is known to be a distance, but the actual distance isn''t known) - ie an arbitrary sized matrix would be one that could be any size.

I think what you are trying to say is that you are making a matrix class for a matrix with a fixed size of 4 x 4 and the values contained by the matrix are not pre-defined.

[btw, I''m the AP who posted earlier]

You are right that in the general case Gaussian elimination is more efficient, however you are using a very small and fixed sized matrix, so I submit that it is likely that the ''classic'' method of inversion is likely to be faster in particular case. This is mostly to do with the reduced number of memory accesses and the lack of need for loops (with the pivot point etc...). I think the minimum number of mults you can get away with in gaussian elimination for a 4 x 4 matrix would be about 50 anyway, and I''m pretty confident I could eliminate 6 multiplies from the ''classic'' method by eliminating common sub-expressions (you''ll find that several of the 2x2 determinates you calculate on the way to the overall 4x4 determinate are repeated several times through out the operation...calc them once & store them in variables, that will probably be more efficient)

Posted 07 March 2004 - 07:33 AM

quote:Original post by optics

then solve for x_{1}to x_{k}. I am having a hard time trying to figure out what b_{1}, b_{2}, b_{3}, b_{4}is suppose to be set to. Do we pull these values out of thin air or is there something I overlooked?

Apply the Gaussain to :

| e11 e12 e13 e14 1 0 0 0 |

| e21 e22 e23 e24 0 1 0 0 |

| e31 e32 e33 e34 0 0 1 0 |

| e41 e42 e43 e44 0 0 0 1 |

Whe it gets reduced to:

| 1 0 0 0 e''11 e''12 e''13 e''14 |

| 0 1 0 0 e''21 e''22 e''23 e''24 |

| 0 0 1 0 e''31 e''32 e''33 e''34 |

| 0 0 0 1 e''41 e''42 e''43 e''44 |

the e'' 4x4 matrix on the right will be your inverse matrix.

Posted 07 March 2004 - 08:58 AM

quote:Original post by optics

then solve for x<sub>1</sub> to x<sub>k</sub>. I am having a hard time trying to figure out what b<sub>1</sub>, b<sub>2</sub>, b<sub>3</sub>, b<sub>4</sub> is suppose to be set to. Do we pull these values out of thin air or is there something I overlooked?quote:

The problem is to solve the system Ax = b

If you calculate the inverse of A, then you can calculate inv(A)b to find x.

You can factor A into upper and lower triangular matrices, to get LUx = b. Then you solve the problem Lc = b, then solve Ux = c. Because the matrices are triangular, the systems can be easily solved with substitution.

You can also attach b as an extra column to A, which is called an augmented matrix, and then use elimination to turn A into an upper triangular matrix, applying the same transformations to b that you do to the rows of A. Then you can solve for x by using the triangular matrix, the new values for b, and substitution.

Posted 07 March 2004 - 11:20 AM

SiCrane: Thank you for the clarification. Now i do have the Gaussian inverse partially implemented I believe but it does not have the identity matrix in the equation. I was wondering two things.

1) Could you show me the mathematical solution to obtain I(A^{-1})?

2)If I post my code would you look over and let me know if I am going in the right direction or am just way off base? I am sure the way that I am doing it is not the most efficient but it seems to be the correct way in my opinion so far.

Then again what do I know?

I thank everyone who has given their help and advice. You guys are really the best.

1) Could you show me the mathematical solution to obtain I(A

2)If I post my code would you look over and let me know if I am going in the right direction or am just way off base? I am sure the way that I am doing it is not the most efficient but it seems to be the correct way in my opinion so far.

Then again what do I know?

I thank everyone who has given their help and advice. You guys are really the best.

Posted 07 March 2004 - 11:27 AM

quote:Original post by optics

1) Could you show me the mathematical solution to obtain I(A^{-1})?

I''m not sure I understand your question. The identity matrix multiplied by any matrix is the other matrix so I(A

quote:

2)If I post my code would you look over and let me know if I am going in the right direction or am just way off base?

Sure.

Posted 07 March 2004 - 11:52 AM

For 1) I was wondering if you might show me the math involved maybe just one step to get the matrix to become what you posted --

Apply the Gaussain to :

Whe it gets reduced to:

2) Here is what I have so far -- be gentle please.

Thanks for the help.

Apply the Gaussain to :

| e11 e12 e13 e14 1 0 0 0 |

| e21 e22 e23 e24 0 1 0 0 |

| e31 e32 e33 e34 0 0 1 0 |

| e41 e42 e43 e44 0 0 0 1 |

Whe it gets reduced to:

| 1 0 0 0 e''11 e''12 e''13 e''14 |

| 0 1 0 0 e''21 e''22 e''23 e''24 |

| 0 0 1 0 e''31 e''32 e''33 e''34 |

| 0 0 0 1 e''41 e''42 e''43 e''44 |

2) Here is what I have so far -- be gentle please.

class CMatrix4x4

{

public:

CMatrix4x4(void); // default constructor

CMatrix4x4(double c11, double c12, double c13, double c14,

double c21, double c22, double c23, double c24,

double c31, double c32, double c33, double c34,

double c41, double c42, double c43, double c44); // constructor

~CMatrix4x4(); // destructor

double Determinant(); // determinant of the matrix

CMatrix4x4& Inverse(); // inverse of the matrix

private:

double e11, e12, e13, e14,

e21, e22, e23, e24,

e31, e32, e33, e34,

e41, e42, e43, e44;

};

double CMatrix4x4::Determinant()

{

return e11*e22*e33*e44 - e11*e22*e34*e43 + e11*e23*e34*e42 - e11*e23*e32*e44

+ e11*e24*e32*e43 - e11*e24*e33*e42 - e12*e23*e34*e41 + e12*e23*e31*e44

- e12*e24*e31*e43 + e12*e24*e33*e41 - e12*e21*e33*e44 + e12*e21*e34*e43

+ e13*e24*e31*e42 - e13*e24*e32*e41 + e13*e21*e32*e44 - e13*e21*e34*e42

+ e13*e22*e34*e41 - e13*e22*e31*e44 - e14*e21*e32*e43 + e14*e21*e33*e42

- e14*e22*e33*e41 + e14*e22*e31*e43 - e14*e23*e31*e42 + e14*e23*e32*e41;

}

CMatrix4x4& CMatrix4x4::Inverse()

{

double s1, s2, s3, s4;

double negX;

if(!this->Determinant()) // if determinant is 0.0 then matrix has no inverse

{ return *this; }

else // otherwise there is an inverse and we need to find it

{

/************************** Step 1 of Gaussian Elimination *******************************/

if( e21 == 1 ) // swap row 2 with row 1 so that e11 == 1

{

s1 = e11; s2 = e12; s3 = e13; s4 = e14;

e11 = e21; e12 = e22; e13 = e23; e14 = e24;

e21 = s1; e22 = s2; e23 = s3; e24 = s4;

}

else if( e31 == 1 ) // swap row 3 with row 1 so that e11 == 1

{

s1 = e11; s2 = e12; s3 = e13; s4 = e14;

e11 = e31; e12 = e32; e13 = e33; e14 = e34;

e31 = s1; e32 = s2; e33 = s3; e34 = s4;

}

else if( e41 == 1 ) // swap row 4 with row 1 so that e11 == 1

{ // do Gaussian step 1

s1 = e11; s2 = e12; s3 = e13; s4 = e14;

e11 = e41; e12 = e42; e13 = e43; e14 = e44;

e41 = s1; e42 = s2; e43 = s3; e44 = s4;

}

// if column 1 does not have a 1 in it then we will do math to make e11 == 1

// sort of hackish but we need to get the code working then we can fix it

//row 2 math

negX = -e21;

e21 = e11/e11*negX + e21;

e22 = e12/e11*negX + e22;

e23 = e13/e11*negX + e23;

e24 = e14/e11*negX + e24;

//row 3 math

negX = -e31;

e31 = e11/e11*negX + e31;

e32 = e12/e11*negX + e32;

e33 = e13/e11*negX + e33;

e34 = e14/e11*negX + e34;

//row 4 math

negX = -e41;

e41 = e11/e11*negX + e41;

e42 = e12/e11*negX + e42;

e43 = e13/e11*negX + e43;

e44 = e14/e11*negX + e44;

/********************** Step 2 of Gaussian Elimination *************************************/

//row 3 math

negX = -e32;

e32 = e22/e22*negX + e32;

e33 = e23/e22*negX + e33;

e34 = e24/e22*negX + e34;

//row 4 math

negX = -e42;

e42 = e22/e22*negX + e42;

e43 = e23/e22*negX + e43;

e44 = e24/e22*negX + e44;

/***************************** Step 3 of Gaussian Elimination *******************************/

//row 4 math

negX = -e43;

e43 = e33/e33*negX + e43;

e44 = e34/e33*negX + e44;

}

return *this;

}

Thanks for the help.

Posted 07 March 2004 - 12:18 PM

For the first part, to get from here:

We can use this (naive) method. First divide the first row by e11.

Then subtract multiples of the first row from each of the second, third and fourth rows.

(e*ij) representing calculations I'm too lazy to type out by hand. Now divide the second row by e*22. Then subtract multiples of row two from frows one, three and four. And so on.

For the second part: In your code, since you're bothering to calculate the determinant, there's no need to do Gaussian reduction - just use the standard determinant method of calculating the inverse. Most of the work is already done if you have the determinant.

edit: formatting

[edited by - SiCrane on March 7, 2004 7:18:40 PM]

to here:

| e11 e12 e13 e14 1 0 0 0 |

| e21 e22 e23 e24 0 1 0 0 |

| e31 e32 e33 e34 0 0 1 0 |

| e41 e42 e43 e44 0 0 0 1 |

| 1 0 0 0 e'11 e'12 e'13 e'14 |

| 0 1 0 0 e'21 e'22 e'23 e'24 |

| 0 0 1 0 e'31 e'32 e'33 e'34 |

| 0 0 0 1 e'41 e'42 e'43 e'44 |

We can use this (naive) method. First divide the first row by e11.

| 1 e12/e11 e13/e11 e14/e11 1/e11 0 0 0 |

| e21 e22 e23 e24 0 1 0 0 |

| e31 e32 e33 e34 0 0 1 0 |

| e41 e42 e43 e44 0 0 0 1 |

Then subtract multiples of the first row from each of the second, third and fourth rows.

| 1 e12/e11 e13/e11 e14/e11 1/e11 0 0 0 |

| 0 e*22 e*23 e*24 -e21/e11 1 0 0 |

| 0 e*32 e*33 e*34 -e31/e11 0 1 0 |

| 0 e*42 e*43 e*44 -e31/e11 0 0 1 |

(e*ij) representing calculations I'm too lazy to type out by hand. Now divide the second row by e*22. Then subtract multiples of row two from frows one, three and four. And so on.

For the second part: In your code, since you're bothering to calculate the determinant, there's no need to do Gaussian reduction - just use the standard determinant method of calculating the inverse. Most of the work is already done if you have the determinant.

edit: formatting

[edited by - SiCrane on March 7, 2004 7:18:40 PM]

Posted 08 March 2004 - 07:39 AM

On the first part, I thank you for that but I also want to ask what you meant by naive method.

Second part: you have to find out of the matrix has a determinant otherwise you will be doing calculations for nothing. If a matrix has a determinant of 0 then there is not an inverse otherwise figure out the inverse. I know this seems backwards to you but I also wanted to code the gaussian elimination method for the knowledge.

I will also end up coding an inverse by way of the determinant method and also by darookie''s suggestion of transposition.

Ok i have done the recode and here is what I have by the math of your suggestion. If its wrong please let me know and I will fix it so that others do not try to copy incorrect code.

Second part: you have to find out of the matrix has a determinant otherwise you will be doing calculations for nothing. If a matrix has a determinant of 0 then there is not an inverse otherwise figure out the inverse. I know this seems backwards to you but I also wanted to code the gaussian elimination method for the knowledge.

I will also end up coding an inverse by way of the determinant method and also by darookie''s suggestion of transposition.

Ok i have done the recode and here is what I have by the math of your suggestion. If its wrong please let me know and I will fix it so that others do not try to copy incorrect code.

class CMatrix4x4

{

public:

CMatrix4x4(void); // default constructor

CMatrix4x4(double c11, double c12, double c13, double c14,

double c21, double c22, double c23, double c24,

double c31, double c32, double c33, double c34,

double c41, double c42, double c43, double c44); // constructor

~CMatrix4x4(); // destructor

Identity(); // identity matrix

double Determinant(); // determinant of the matrix

CMatrix4x4& Inverse(); // inverse of the matrix

private:

double e11, e12, e13, e14,

e21, e22, e23, e24,

e31, e32, e33, e34,

e41, e42, e43, e44;

};

CMatrix4x4::Identity()

{

e11 = e22 = e33 = e44 = 1.0f;

e12 = e13 = e14 = e21 = e23 = e24 =

e31 = e32 = e34 = e41 = e42 = e43 = 0.0f;

}

double CMatrix4x4::Determinant()

{

return e11*e22*e33*e44 - e11*e22*e34*e43 + e11*e23*e34*e42 - e11*e23*e32*e44

+ e11*e24*e32*e43 - e11*e24*e33*e42 - e12*e23*e34*e41 + e12*e23*e31*e44

- e12*e24*e31*e43 + e12*e24*e33*e41 - e12*e21*e33*e44 + e12*e21*e34*e43

+ e13*e24*e31*e42 - e13*e24*e32*e41 + e13*e21*e32*e44 - e13*e21*e34*e42

+ e13*e22*e34*e41 - e13*e22*e31*e44 - e14*e21*e32*e43 + e14*e21*e33*e42

- e14*e22*e33*e41 + e14*e22*e31*e43 - e14*e23*e31*e42 + e14*e23*e32*e41;

}

CMatrix4x4& CMatrix4x4::Inverse()

{

double s1, s2, s3, s4;

double negX, div;

CMatrix4x4 m4x4Inv;

m4x4Inv.Identity();

if(!this->Determinant()) // if determinant is 0.0 then matrix has no inverse

{ return *this; }

else // otherwise there is an inverse and we need to find it

{

/************************** Step 1 of Gaussian Elimination *******************************/

if( e21 == 1 ) // swap row 2 with row 1 so that e11 == 1

{

s1 = e11; s2 = e12; s3 = e13; s4 = e14;

e11 = e21; e12 = e22; e13 = e23; e14 = e24;

e21 = s1; e22 = s2; e23 = s3; e24 = s4;

}

else if( e31 == 1 ) // swap row 3 with row 1 so that e11 == 1

{

s1 = e11; s2 = e12; s3 = e13; s4 = e14;

e11 = e31; e12 = e32; e13 = e33; e14 = e34;

e31 = s1; e32 = s2; e33 = s3; e34 = s4;

}

else if( e41 == 1 ) // swap row 4 with row 1 so that e11 == 1

{

s1 = e11; s2 = e12; s3 = e13; s4 = e14;

e11 = e41; e12 = e42; e13 = e43; e14 = e44;

e41 = s1; e42 = s2; e43 = s3; e44 = s4;

}

//if column 1 does not have a 1 in it then we will do math to make e11 == 1

//sort of hackish but we need to get the code working then we can fix it

// row 1 math

div = e11;

e11 = e11/div;

e12 = e12/div;

e13 = e13/div;

e14 = e14/div;

m4x4Inv.e11 = 1/div;

// row 2 math

negX = -e21;

e21 = e11*negX + e21;

e22 = e12*negX + e22;

e23 = e13*negX + e23;

e24 = e14*negX + e24;

m4x4Inv.e21 = negX/div;

// row 3 math

negX = -e31;

e31 = negX + e31;

e32 = e12*negX + e32;

e33 = e13*negX + e33;

e34 = e14*negX + e34;

m4x4Inv.e31 = negX/div;

// row 4 math

negX = -e41;

e41 = negX + e41;

e42 = e12*negX + e42;

e43 = e13*negX + e43;

e44 = e14*negX + e44;

m4x4Inv.e41 = negX/div;

/********************* Step 2 of Gaussian Elimination *************************************/

div = e22;

// row 2 math

e22 = e22/div;

e23 = e23/div;

e24 = e24/div;

m4x4Inv.e22 = 1/div;

// row 1 math

negX = -e12;

e12 = negX + e12;

e13 = e23*negX + e13;

e14 = e24*negX + e14;

m4x4Inv.e12 = negX/div;

// row 3 math

negX = -e32;

e32 = negX + e32;

e33 = e23*negX + e33;

e34 = e24*negX + e34;

m4x4Inv.e32 = negX/div;

// row 4 math

negX = -e42;

e42 = negX + e42;

e43 = e23*negX + e43;

e44 = e24*negX + e44;

m4x4Inv.e42 = negX/div;

/***************************** Step 3 of Gaussian Elimination *******************************/

// row 3 math

div = e33;

e33 = e33/div;

e34 = e34/div;

m4x4Inv.e33 = 1/div;

//row 1 math

negX = -e13;

e13 = negX + e13;

e14 = e34*negX + e14;

m4x4Inv.e13 = negX/div;

// row 2 math

negX = -e23;

e23 = negX + e23;

e24 = e34*negX + e24;

m4x4Inv.e23 = negX/div;

// row 4 math

negX = -e43;

e43 = negX + e43;

e44 = e34*negX + e44;

m4x4Inv.e43 = negX/div;

/***************************** Step 4 of Gaussian Elimination *******************************/

// row 4 math

div = e44;

e44 = e44/div;

m4x4Inv.e44 = 1/div;

// row 1 math

negX = -e14;

e14 = negX + e14;

m4x4Inv.e14 = negX/div;

// row 2 math

negX = -e24;

e24 = negX + e24;

m4x4Inv.e24 = negX/div;

// row 3 math

negX = -e34;

e34 = negX + e34;

m4x4Inv.e34 = negX/div;

// assign m4x4Inv to the matrix element by element

e11 = m4x4Inv.e11; e12 = m4x4Inv.e12; e13 = m4x4Inv.e13; e14 = m4x4Inv.e14;

e21 = m4x4Inv.e21; e22 = m4x4Inv.e22; e23 = m4x4Inv.e23; e24 = m4x4Inv.e24;

e31 = m4x4Inv.e31; e32 = m4x4Inv.e32; e33 = m4x4Inv.e33; e34 = m4x4Inv.e34;

e41 = m4x4Inv.e41; e42 = m4x4Inv.e42; e43 = m4x4Inv.e43; e44 = m4x4Inv.e44;

}

return *this;

}

Posted 08 March 2004 - 09:12 AM

quote:Original post by optics

On the first part, I thank you for that but I also want to ask what you meant by naive method.

For one thing, I didn''t bother to check if e11 was 0 or not before I divided the row by it. Nor did I check to see if e*22 was 0 or not before I divided row 2 by it. Also because of the nature of floating point arithmetic, it''s better, instead of just choosing the number in the first row to divide by, swap the rows so that you divide by the number with the largest value in the first column, etc.

quote:

Second part: you have to find out of the matrix has a determinant otherwise you will be doing calculations for nothing.

Well, since you''re throwing away the determinant calculation, you''re already doing calculations for nothing. If you''re doing Gaussian elimination, then you should just procede with the elimination steps until you get to a point where the you don''t have a valid divisor.

As for the actual code, there are a couple issues with the algorithm implementation.

1: you aren''t checking for divide by zero.

2: when you do the computations for the second, third and fourth rows, you don''t update the previous columns in the inverse half of the augmented matrix.

Also purely, in terms of C++ coding style there are some additional issues.

1: your code would look a lot prettier if you you used std::swap() instead of doing manual swaps.

2: the code would also be a lot cleaner if you create a row vector class to encapsulate some of the matrix math going on in the elimination routine.

3: Some of the function names are misleading. Like the Identity() function. It turns the matrix into an identity matrix, but the name sounds like it might return an identity matrix. Similarly the Inverse() function sounds like it should return the matrix inverse instead of making the current matrix into its inverse.

Posted 08 March 2004 - 11:31 AM

SiCrane: Thank you for all the tips.

Todo for the algorigthm: I know I am not checking for divide by zero and I appreciate you pointing that out to me because knowing me I would forget it in the end. But I also plan on posting the finished code for more critiques.

I also thought it was better if e11 == 1 rather than e11 being the largest value in the column. I have not studied this evaluation at all so I would not know best. If you have a link about this I would be more than willing to read and learn from it. It would also be most appreciated.

Could you show the math for step 2 on updating the previous columns? I have tried figuring this out but it looks to me like I am just multiplying the previous columns by -1.

for example for row 3 in step 2

should I also be adding m4x4Inv.e_{ij} in this case m4x4Inv.e31 to the same element? Or am I just way off.

I know it might be asking a lot of you but could you write out the math of that step if only for one update so that I can see it?

Also in the case of a division by zero should I just set those elements to zero?

Thank you for the help on this so far.

Todo for the algorigthm: I know I am not checking for divide by zero and I appreciate you pointing that out to me because knowing me I would forget it in the end. But I also plan on posting the finished code for more critiques.

I also thought it was better if e11 == 1 rather than e11 being the largest value in the column. I have not studied this evaluation at all so I would not know best. If you have a link about this I would be more than willing to read and learn from it. It would also be most appreciated.

Could you show the math for step 2 on updating the previous columns? I have tried figuring this out but it looks to me like I am just multiplying the previous columns by -1.

for example for row 3 in step 2

/** Step 2 of Gaussian Elimination **/

// row 3 math

// div = e22

negX = -e32;

e32 = negX + e32;

e33 = e23*negX + e33;

e34 = e24*negX + e34;

m4x4Inv.e31 = negX/div*m4x4Inv.e31;

m4x4Inv.e32 = negX/div;

should I also be adding m4x4Inv.e

I know it might be asking a lot of you but could you write out the math of that step if only for one update so that I can see it?

Also in the case of a division by zero should I just set those elements to zero?

Thank you for the help on this so far.