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.

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

D3DXMATRIX4

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

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.

I just want to thank you all for helping me sort this problem out and I hope that others will learn from the code that I have written.

SiCrane: If you are ever coming to Dallas, Texas let me know and I will have to buy you a case of beer or three for your extensive help on this little project. Without your help I don''t think I could this done.

joanusdmentia: If your algorithm is quicker and/or better than the one I have put here would you mind sharing it with us? We are always in a constant state of learning and if your implementation provides a faster/better way of doing this then I for one would love to see it.

psamty10: I was looking more on information of inverting a matrix by Gaussian Elimination but thank for the link maybe someone else can use that companies code.

I will be coding a Vector Class that will replace the double elements in this Matrix Class and post the code either in this thread or I will start a new one.
The final step will be to make this and the Vector Class templated for data types and Size( i. e. NxN )

Thank you all again and I hope to maybe see my code in some working games with proper credit of course or if modified a simple hey would be enough for me.

Thanks,

Logan Watt

I realize I''m getting to this thread late, but I just wanted to toss in my two cents on 4x4 matrix inversion.

I use a Cramer''s Rule inverter. You can see it here: http://aoeu.snth.net/vec/vec.c

The "m_inv" function performs a matrix inverse in a completely unrolled, nicely pipelineable manner. It performs well. The code uses an array of doubles rather than psuedo-subscript variable names, but you can fix that with some global replaces. The double array is used because this is the array format returned by OpenGL matrix queries. You might also want to add a zero test.