• ### Announcements

#### Archived

This topic is now archived and is closed to further replies.

# Invert 4x4 Matrix Gaussian Elimination

## Recommended Posts

optics    122
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]

##### Share on other sites
SiCrane    11839
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.

##### Share on other sites
Manbot1    122
We have to solve a 9x9 matrix by next Friday for class, I''ll tell you how it goes.

##### Share on other sites
optics    122
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.

##### Share on other sites
SiCrane    11839
If you end with a row of all 0s then the matrix is non-invertible. So that really isn''t an issue.

##### Share on other sites
darookie    1441
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.

##### Share on other sites
Guest Anonymous Poster
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:
(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...

##### Share on other sites
optics    122
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]

##### Share on other sites
joanusdmentia    1060
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]

##### Share on other sites
optics    122
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 b1 ... bk 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 x1 to xk. I am having a hard time trying to figure out what b1, b2, b3, b4 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]

##### Share on other sites
Guest Anonymous Poster
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)

##### Share on other sites
SiCrane    11839
quote:
Original post by optics

then solve for x1 to xk. I am having a hard time trying to figure out what b1, b2, b3, b4 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.

##### Share on other sites
Guest Anonymous Poster
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.

##### Share on other sites
optics    122
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.

##### Share on other sites
SiCrane    11839
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-1) would just be A-1.
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.

##### Share on other sites
optics    122
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 :
      | 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 matrixprivate:	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.

##### Share on other sites
SiCrane    11839
For the first part, to get from 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 |
to here:
      | 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]

##### Share on other sites
optics    122
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.

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

##### Share on other sites
SiCrane    11839
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.

##### Share on other sites
optics    122
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

/** 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.eij 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.

##### Share on other sites
SiCrane    11839
Let me demonstrate with numerical examples, since writing out the symbols would be too much of a pain. I''m going to dump this in a source box, because step by step a Gaussian elimination takes up a lot of room.
Here''s the beginning augmented matrix. The matrix to be invertedis on the left, and the identity matrix is on the right.-2.000	0.000	3.000	4.000	1.000	0.000	0.000	0.000	1.000	0.000	1.000	-2.000	0.000	1.000	0.000	0.000	-1.000	5.000	4.000	5.000	0.000	0.000	1.000	0.000	-2.000	3.000	8.000	7.000	0.000	0.000	0.000	1.000	Step 1: First we find the pivot, which we determine by lookingat the absolute value of the first each element in the firstcolumn. In this case we don''t do anything because the firstelement in the first column has an absolute value greater thanor equal to all the other elements in the first column.-2.000	0.000	3.000	4.000	1.000	0.000	0.000	0.000	1.000	0.000	1.000	-2.000	0.000	1.000	0.000	0.000	-1.000	5.000	4.000	5.000	0.000	0.000	1.000	0.000	-2.000	3.000	8.000	7.000	0.000	0.000	0.000	1.000	Then we divide every element in the first row by the value ofthe first element in the first row.1.000	0.000	-1.500	-2.000	-0.500	0.000	0.000	0.000	1.000	0.000	1.000	-2.000	0.000	1.000	0.000	0.000	-1.000	5.000	4.000	5.000	0.000	0.000	1.000	0.000	-2.000	3.000	8.000	7.000	0.000	0.000	0.000	1.000	Now we subtract 1 times row 1 from row 2.1.000	0.000	-1.500	-2.000	-0.500	0.000	0.000	0.000	0.000	0.000	2.500	0.000	0.500	1.000	0.000	0.000	-1.000	5.000	4.000	5.000	0.000	0.000	1.000	0.000	-2.000	3.000	8.000	7.000	0.000	0.000	0.000	1.000	Now we subtract -1 times row from row 3.1.000	0.000	-1.500	-2.000	-0.500	0.000	0.000	0.000	0.000	0.000	2.500	0.000	0.500	1.000	0.000	0.000	0.000	5.000	2.500	3.000	-0.500	0.000	1.000	0.000	-2.000	3.000	8.000	7.000	0.000	0.000	0.000	1.000	Now we subtract -2 times row 1 from row 4.1.000	0.000	-1.500	-2.000	-0.500	0.000	0.000	0.000	0.000	0.000	2.500	0.000	0.500	1.000	0.000	0.000	0.000	5.000	2.500	3.000	-0.500	0.000	1.000	0.000	0.000	3.000	5.000	3.000	-1.000	0.000	0.000	1.000	Step 2: Again find the pivot, in this case we''re looking atcolumn 2. We can''t use the 0 value in row two, so we need to swaprow 2 with another row. Row 3 has the largest value for column 2so we swap rows 2 and 3.1.000	0.000	-1.500	-2.000	-0.500	0.000	0.000	0.000	0.000	5.000	2.500	3.000	-0.500	0.000	1.000	0.000	0.000	0.000	2.500	0.000	0.500	1.000	0.000	0.000	0.000	3.000	5.000	3.000	-1.000	0.000	0.000	1.000	Now divide row 3 by the pivot value1.000	0.000	-1.500	-2.000	-0.500	0.000	0.000	0.000	0.000	1.000	0.500	0.600	-0.100	0.000	0.200	0.000	0.000	0.000	2.500	0.000	0.500	1.000	0.000	0.000	0.000	3.000	5.000	3.000	-1.000	0.000	0.000	1.000	Since row 1 has 0 in column 2, we don''t need to do anything torow 1.1.000	0.000	-1.500	-2.000	-0.500	0.000	0.000	0.000	0.000	1.000	0.500	0.600	-0.100	0.000	0.200	0.000	0.000	0.000	2.500	0.000	0.500	1.000	0.000	0.000	0.000	3.000	5.000	3.000	-1.000	0.000	0.000	1.000	Since row 3 has 0 in column 2, we don''t need to do anything torow 3.1.000	0.000	-1.500	-2.000	-0.500	0.000	0.000	0.000	0.000	1.000	0.500	0.600	-0.100	0.000	0.200	0.000	0.000	0.000	2.500	0.000	0.500	1.000	0.000	0.000	0.000	3.000	5.000	3.000	-1.000	0.000	0.000	1.000	Subtract 3 times row 2 from row 4.1.000	0.000	-1.500	-2.000	-0.500	0.000	0.000	0.000	0.000	1.000	0.500	0.600	-0.100	0.000	0.200	0.000	0.000	0.000	2.500	0.000	0.500	1.000	0.000	0.000	0.000	0.000	3.500	1.200	-0.700	0.000	-0.600	1.000	Step 3: No we''re working on column 3. Row 4 has a larger value incolumn 3 than row 3, so we swap rows 3 and 4.1.000	0.000	-1.500	-2.000	-0.500	0.000	0.000	0.000	0.000	1.000	0.500	0.600	-0.100	0.000	0.200	0.000	0.000	0.000	3.500	1.200	-0.700	0.000	-0.600	1.000	0.000	0.000	2.500	0.000	0.500	1.000	0.000	0.000	Then divide row 3 by it''s value in column 3.1.000	0.000	-1.500	-2.000	-0.500	0.000	0.000	0.000	0.000	1.000	0.500	0.600	-0.100	0.000	0.200	0.000	0.000	0.000	1.000	0.343	-0.200	0.000	-0.171	0.286	0.000	0.000	2.500	0.000	0.500	1.000	0.000	0.000	Subtract -1.5 times row 3 from row 1.1.000	0.000	0.000	-1.486	-0.800	0.000	-0.257	0.429	0.000	1.000	0.500	0.600	-0.100	0.000	0.200	0.000	0.000	0.000	1.000	0.343	-0.200	0.000	-0.171	0.286	0.000	0.000	2.500	0.000	0.500	1.000	0.000	0.000	Then subtract 0.5 times row 3 from row 2.1.000	0.000	0.000	-1.486	-0.800	0.000	-0.257	0.429	0.000	1.000	0.000	0.429	-0.000	0.000	0.286	-0.143	0.000	0.000	1.000	0.343	-0.200	0.000	-0.171	0.286	0.000	0.000	2.500	0.000	0.500	1.000	0.000	0.000	Then subtract 2.5 times row 3 from row 41.000	0.000	0.000	-1.486	-0.800	0.000	-0.257	0.429	0.000	1.000	0.000	0.429	-0.000	0.000	0.286	-0.143	0.000	0.000	1.000	0.343	-0.200	0.000	-0.171	0.286	0.000	0.000	0.000	-0.857	1.000	1.000	0.429	-0.714	Step 4: Even though row 1 has a greater absolute value in thefourth column than row 4, we don''t do any swapping because row1 has already been processed. In fact, this is the last row sono swapping can be done. So we just divide row 4 by -0.8571.000	0.000	0.000	-1.486	-0.800	0.000	-0.257	0.429	0.000	1.000	0.000	0.429	-0.000	0.000	0.286	-0.143	0.000	0.000	1.000	0.343	-0.200	0.000	-0.171	0.286	0.000	0.000	0.000	1.000	-1.167	-1.167	-0.500	0.833	Then subtract -1.486 times row 4 from row 1.1.000	0.000	0.000	0.000	-2.533	-1.733	-1.000	1.667	0.000	1.000	0.000	0.429	-0.000	0.000	0.286	-0.143	0.000	0.000	1.000	0.343	-0.200	0.000	-0.171	0.286	0.000	0.000	0.000	1.000	-1.167	-1.167	-0.500	0.833	Then subtract .429 times row 4 from row 2.1.000	0.000	0.000	0.000	-2.533	-1.733	-1.000	1.667	0.000	1.000	0.000	0.000	0.500	0.500	0.500	-0.500	0.000	0.000	1.000	0.343	-0.200	0.000	-0.171	0.286	0.000	0.000	0.000	1.000	-1.167	-1.167	-0.500	0.833	Then subtract .343 times row 4 from row 3.1.000	0.000	0.000	0.000	-2.533	-1.733	-1.000	1.667	0.000	1.000	0.000	0.000	0.500	0.500	0.500	-0.500	0.000	0.000	1.000	0.000	0.200	0.400	0.000	0.000	0.000	0.000	0.000	1.000	-1.167	-1.167	-0.500	0.833	Now the matrix has been reduced so the 4x4 on the right is(roughly) equivalent to the inverse of the original matrix.

If at any point all the unprocessed rows had 0''s for elements in the current column being worked on, then the matrix would have been non-invertible, and your function should react accordingly.

##### Share on other sites
optics    122
SiCrane minus the check for division by zero and the check on each step for eij == 0 where i == j, this should be completed.
I should be putting error checking in shortly.

So that takes care of the 1) problems you had with the algorithm.

2) For the C++ I changed the name of the Identity function to make it more to standards I guess you could say. Any other suggestions you could give on this subject?

I do plan on coding a Vector Class and make those part of the matrix class by inheritance or such. I should probably remove the other code as it is wrong or if anyone posts to keep it here I will.

Here is the Finished code w/o error checking.
#include <iostream>#include <cmath> // for fabs();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	MakeIdentity(); // turn matrix into identity matrix	CMatrix4x4& Inverse(); // inverse of the matrixprivate:	double e11, e12, e13, e14, 		  e21, e22, e23, e24,		  e31, e32, e33, e34,		  e41, e42, e43, e44;};CMatrix4x4::MakeIdentity(){	e11 = e22 = e33 = e44 = 1.0f;			e12 = e13 = e14 = e21 = e23 = e24 = 	e31 = e32 = e34 = e41 = e42 = e43 = 0.0f;}CMatrix4x4& CMatrix4x4::Inverse(){	double s1, s2, s3, s4;	double negX, div;	double r1fab, r2fab, r3fab, r4fab;	CMatrix4x4 m4x4Inv;	m4x4Inv.MakeIdentity();/*========================= Step 1 of Gaussian Elimination ==============================*/		r1fab = fabs(e11);		r2fab = fabs(e21);		r3fab = fabs(e31);		r4fab = fabs(e41);		// Pivot Point Test on Column 1				if( r2fab > r1fab ) 		{						s1 = e11;  			s2 = e12;  			s3 = e13;  			s4 = e14;			e11 = e21; 			e12 = e22; 			e13 = e23; 			e14 = e24;			e21 = s1;  			e22 = s2;  			e23 = s3;  			e24 = s4;			 			s1 = m4x4Inv.e11; 			s2 = m4x4Inv.e12; 			s3 = m4x4Inv.e13; 			s4 = m4x4Inv.e14;			m4x4Inv.e11 = m4x4Inv.e21; 			m4x4Inv.e12 = m4x4Inv.e22; 			m4x4Inv.e13 = m4x4Inv.e23;			m4x4Inv.e14 = m4x4Inv.e24;			m4x4Inv.e21 = s1;			m4x4Inv.e22 = s2;			m4x4Inv.e23 = s3;			m4x4Inv.e24 = s4;		}		if( r3fab > r1fab ) 		{				s1 = e11; 			s2 = e12; 			s3 = e13; 			s4 = e14;			e11 = e31; 			e12 = e32; 			e13 = e33; 			e14 = e34;			e31 = s1; 			e32 = s2; 			e33 = s3; 			e34 = s4;			s1 = m4x4Inv.e11; 			s2 = m4x4Inv.e12; 			s3 = m4x4Inv.e13; 			s4 = m4x4Inv.e14;			m4x4Inv.e11 = m4x4Inv.e31; 			m4x4Inv.e12 = m4x4Inv.e32; 			m4x4Inv.e13 = m4x4Inv.e33;			m4x4Inv.e14 = m4x4Inv.e34;			m4x4Inv.e31 = s1;			m4x4Inv.e32 = s2;			m4x4Inv.e33 = s3;			m4x4Inv.e34 = s4;		}		if( r4fab > r1fab ) 		{								s1 = e11; 			s2 = e12; 			s3 = e13; 			s4 = e14;			e11 = e41; 			e12 = e42;			e13 = e43; 			e14 = e44;			e41 = s1; 			e42 = s2; 			e43 = s3; 			e44 = s4;			s1 = m4x4Inv.e11; 			s2 = m4x4Inv.e12; 			s3 = m4x4Inv.e13; 			s4 = m4x4Inv.e14;			m4x4Inv.e11 = m4x4Inv.e41; 			m4x4Inv.e12 = m4x4Inv.e42; 			m4x4Inv.e13 = m4x4Inv.e43;			m4x4Inv.e14 = m4x4Inv.e44;			m4x4Inv.e41 = s1;			m4x4Inv.e42 = s2;			m4x4Inv.e43 = s3;			m4x4Inv.e44 = s4;		}				/*if column 1 was not pivoted then do Step 1 in Gaussian Elimination 	        			 row 1 math			 Divide first row by e11 to make e11 == 1			 */				div = e11;				e11 /= div;				e12 /= div;				e13 /= div;				e14 /= div;				m4x4Inv.e11 = m4x4Inv.e11/div;				m4x4Inv.e12 = m4x4Inv.e12/div;				m4x4Inv.e13 = m4x4Inv.e13/div;				m4x4Inv.e14 = m4x4Inv.e14/div;				// row 2 math					negX = -e21;					e21 += e11*negX;					e22 += e12*negX;					e23 += e13*negX;					e24 += e14*negX;					m4x4Inv.e21 += m4x4Inv.e11*negX;					m4x4Inv.e22 += m4x4Inv.e12*negX;					m4x4Inv.e23 += m4x4Inv.e13*negX;					m4x4Inv.e24 += m4x4Inv.e14*negX;				// row 3 math					negX = -e31;					e31 += negX;					e32 += e12*negX;					e33 += e13*negX;					e34 += e14*negX;					m4x4Inv.e31 += m4x4Inv.e11*negX;					m4x4Inv.e32 += m4x4Inv.e12*negX;					m4x4Inv.e33 += m4x4Inv.e13*negX;					m4x4Inv.e34 += m4x4Inv.e14*negX;							// row 4 math					negX = -e41;					e41 += negX;					e42 += e12*negX;					e43 += e13*negX;					e44 += e14*negX;					m4x4Inv.e41 += m4x4Inv.e11*negX;					m4x4Inv.e42 += m4x4Inv.e12*negX;					m4x4Inv.e43 += m4x4Inv.e13*negX;					m4x4Inv.e44 += m4x4Inv.e14*negX;/*======================== Step 2 of Gaussian Elimination ================================*/			// Pivot Point Test for Column 2 excluding row 1		r2fab = fabs(e22);		r3fab = fabs(e32);		r4fab = fabs(e42);		if( r3fab > r2fab ) 		{				s1 = e21; 			s2 = e22; 			s3 = e23; 			s4 = e24;			e21 = e31; 			e22 = e32; 			e23 = e33; 			e24 = e34;			e31 = s1; 			e32 = s2; 			e33 = s3; 			e34 = s4;			s1 = m4x4Inv.e21; 			s2 = m4x4Inv.e22; 			s3 = m4x4Inv.e23; 			s4 = m4x4Inv.e24;			m4x4Inv.e21 = m4x4Inv.e31; 			m4x4Inv.e22 = m4x4Inv.e32; 			m4x4Inv.e23 = m4x4Inv.e33;			m4x4Inv.e24 = m4x4Inv.e34;			m4x4Inv.e31 = s1;			m4x4Inv.e32 = s2;			m4x4Inv.e33 = s3;			m4x4Inv.e34 = s4;		}		if( r4fab > r2fab ) 		{								s1 = e21; 			s2 = e22; 			s3 = e23; 			s4 = e24;			e21 = e41; 			e22 = e42;			e23 = e43; 			e24 = e44;			e41 = s1; 			e42 = s2; 			e43 = s3; 			e44 = s4;			s1 = m4x4Inv.e21; 			s2 = m4x4Inv.e22; 			s3 = m4x4Inv.e23; 			s4 = m4x4Inv.e24;			m4x4Inv.e21 = m4x4Inv.e41; 			m4x4Inv.e22 = m4x4Inv.e42; 			m4x4Inv.e23 = m4x4Inv.e43;			m4x4Inv.e24 = m4x4Inv.e44;			m4x4Inv.e41 = s1;			m4x4Inv.e42 = s2;			m4x4Inv.e43 = s3;			m4x4Inv.e44 = s4;		}						/* row 2 math			   Divide row 2 by e22 to make e22 == 1			*/				div = e22;				e22 /= div;				e23 /= div;				e24 /= div;				m4x4Inv.e21 /= div;				m4x4Inv.e22 /= div;				m4x4Inv.e23 /= div;				m4x4Inv.e24 /= div;				// row 1 math					negX = -e12;					e12 += negX;					e13 += e23*negX;					e14 += e24*negX;					m4x4Inv.e11 += negX*m4x4Inv.e21;					m4x4Inv.e12 += negX*m4x4Inv.e22;					m4x4Inv.e13 += negX*m4x4Inv.e23;					m4x4Inv.e14 += negX*m4x4Inv.e24;							// row 3 math					negX = -e32;					e32 += negX;					e33 += e23*negX;					e34 += e24*negX;					m4x4Inv.e31 += negX*m4x4Inv.e21;					m4x4Inv.e32 += negX*m4x4Inv.e22;					m4x4Inv.e33 += negX*m4x4Inv.e23;					m4x4Inv.e34 += negX*m4x4Inv.e24;							// row 4 math					negX = -e42;					e42 += negX;					e43 += e23*negX;					e44 += e24*negX;					m4x4Inv.e41 += negX*m4x4Inv.e21;					m4x4Inv.e42 += negX*m4x4Inv.e22;					m4x4Inv.e43 += negX*m4x4Inv.e23;					m4x4Inv.e44 += negX*m4x4Inv.e24;/*====================== Step 3 of Gaussian Elimination ================================*/			r3fab = fabs(e33);			r4fab = fabs(e43);		// Pivot Point Test for Column 3 exluding row 1 and 2		if( r4fab > r3fab ) 		{								s1 = e31; 			s2 = e32; 			s3 = e33; 			s4 = e34;			e31 = e41; 			e32 = e42;			e33 = e43; 			e34 = e44;			e41 = s1; 			e42 = s2; 			e43 = s3; 			e44 = s4;			s1 = m4x4Inv.e31; 			s2 = m4x4Inv.e32; 			s3 = m4x4Inv.e33; 			s4 = m4x4Inv.e34;			m4x4Inv.e31 = m4x4Inv.e41; 			m4x4Inv.e32 = m4x4Inv.e42; 			m4x4Inv.e33 = m4x4Inv.e43;			m4x4Inv.e34 = m4x4Inv.e44;			m4x4Inv.e41 = s1;			m4x4Inv.e42 = s2;			m4x4Inv.e43 = s3;			m4x4Inv.e44 = s4;		}			// row 3 math			div = e33;			e33 /= div;			e34 /= div;			m4x4Inv.e31 /= div;			m4x4Inv.e32 /= div;			m4x4Inv.e33 /= div;			m4x4Inv.e34 /= div;			//row 1 math				negX = -e13;				e13 += negX;				e14 += e34*negX;				m4x4Inv.e11 += negX*m4x4Inv.e31;				m4x4Inv.e12 += negX*m4x4Inv.e32;				m4x4Inv.e13 += negX*m4x4Inv.e33;				m4x4Inv.e14 += negX*m4x4Inv.e34;						// row 2 math				negX = -e23;				e23 += negX;				e24 += e34*negX;				m4x4Inv.e21 += negX*m4x4Inv.e31;				m4x4Inv.e22 += negX*m4x4Inv.e32;				m4x4Inv.e23 += negX*m4x4Inv.e33;				m4x4Inv.e24 += negX*m4x4Inv.e34;			// row 4 math				negX = -e43;				e43 += negX;				e44 += e34*negX;				m4x4Inv.e41 += negX*m4x4Inv.e31;				m4x4Inv.e42 += negX*m4x4Inv.e32;				m4x4Inv.e43 += negX*m4x4Inv.e33;				m4x4Inv.e44 += negX*m4x4Inv.e34;/*======================= Step 4 of Gaussian Elimination ===============================*/						// No Pivot Point Test since we are at row 4			// row 4 math			div = e44;			e44 /= div;			m4x4Inv.e41 /= div;			m4x4Inv.e42 /= div;			m4x4Inv.e43 /= div;			m4x4Inv.e44 /= div;			// row 1 math				negX = -e14;				e14 += negX;				m4x4Inv.e11 += negX*m4x4Inv.e41;				m4x4Inv.e12 += negX*m4x4Inv.e42;				m4x4Inv.e13 += negX*m4x4Inv.e43;				m4x4Inv.e14 += negX*m4x4Inv.e44;			// row 2 math				negX = -e24;				e24 = negX + e24;				m4x4Inv.e21 += negX*m4x4Inv.e41;				m4x4Inv.e22 += negX*m4x4Inv.e42;				m4x4Inv.e23 += negX*m4x4Inv.e43;				m4x4Inv.e24 += negX*m4x4Inv.e44;			// row 3 math				negX = -e34;				e34 = negX + e34;				m4x4Inv.e31 += negX*m4x4Inv.e41;				m4x4Inv.e32 += negX*m4x4Inv.e42;				m4x4Inv.e33 += negX*m4x4Inv.e43;				m4x4Inv.e34 += negX*m4x4Inv.e44;			std::swap(*this, m4x4Inv);			  return *this;}

psamty10    148
D3DXMATRIX4

##### Share on other sites
psamty10    148
since my first reply was so stupid, i thought id make my stupid self useful.
Magic Software has a good implementation at Clicky

Perhaps that will give you ze idea

##### Share on other sites
SiCrane    11839
I looked the code over and ran a few test cases on it, and I don''t see any problems from a program correctness point of view (other than the ones that you already mentioned.)

From a C++ coding style point of view, there are a few issues. At least in my opinion, Inverse() still isn''t a particularly good name. From the name alone it sounds like the function should return the inverse of the matrix. Instead it inverts the matrix. As a general rule of thumb, functions with noun-names shouldn''t directly modify the class they are members of, especially if they don''t have any function arguments.

Also when you swap rows, you swap individual elements using temporary variables. This could be changed to using std::swap, which will reduce the length of the function by quite a bit.

Again, using row vectors instead of storing the matrix elements individually would also probably reduce the code size even more. At least it would make the inverse function more managable.

And this is a minor thing: some of the indentation seems a little inconsistent. Like for the fabs() assignment statements.