Matrix inverse

[ L  T ] -1    [ L'  -L'*T ][      ]     = [           ]   with L' := L^-1[ 0  1 ]       [  0    1   ]with the proove[ L  T ]   [ L'  -L'*T ]   [ L*L'  L*(-L'*T)+T ]   [ I  -(L*L')*T+T ]   [ I  0 ][      ] * [           ] = [                   ] = [                ] = [      ]  q.e.d.[ 0  1 ]   [  0    1   ]   [  0          1     ]   [ 0        1     ]   [ 0  1 ]

result._rotating = _rotating;result._scaling = _scaling;if(_rotating && _scaling) {	// presetting linear sub-matrix w/ adjoint of linear sub-matrix of self	result._b0 = _b1.cross(_b2);	result._b1 = _b2.cross(_b0);	result._b2 = _b0.cross(_b1);	// scaling adjoint sub-matrix by inverse of determinant results in inverse sub-matrix; the determinant is expected	// to never be zero since that would not represent a valid affine transformation	const register scalar_g inverseOfDet = 1.0f/_b0.dot(result._b0);	result._b0.multiply(inverseOfDet);	result._b1.multiply(inverseOfDet);	result._b2.multiply(inverseOfDet);} else if(_rotating) {	// rotating only: inverse is identical to transpose	result._b0.set(_b0.s0,_b1.s0,_b2.s0);	result._b1.set(_b0.s1,_b1.s1,_b2.s1);	result._b2.set(_b0.s2,_b1.s2,_b2.s2);} else if(_scaling) {	// scaling only: matrix is diagonal, and so is its inverse; assuming that there is no zero	result._b0.set(scalar_g(1)/_b0.s0,0,0);	result._b1.set(0,scalar_g(1)/_b1.s1,0);	result._b2.set(0,0,scalar_g(1)/_b2.s2);} else {	// identity	result._b0 = Vector_t::UNIT_0;	result._b1 = Vector_t::UNIT_1;	result._b2 = Vector_t::UNIT_2;}

/** * This function sets 'this' matrix to the inverse of 'this' matrix. * \return 'true' if the matrix has an inverse and it was set to the inverse, else 'false' * and the matrix isn't changed. */bool CMatrix::setInverse(){	// first split the matrix using the top row and each column to find the determinant	// save the values for the inverse	float temp[16];		// temp for the inverse matrix	float deter = 0.0f;	// determinant for the matrix	// calculate the determinant of the matrix (top row of the matrix)	float m[9];	// temp 3x3 submatrix	// keep in mind deter * (-1)^(j+1) where j is the column that is being ignored	split(m, 0, 0);	temp[0] = -calcDeter(m, 3);	// first deter bit * (-1)^1	split(m, 1, 0);	temp[1] = calcDeter(m, 3);	// second deter bit * (-1)^(1+1)			split(m, 2, 0);	temp[2] = -calcDeter(m, 3);	// third deter bit * (-1)^(2+1)	split(m, 3, 0);	temp[3] = calcDeter(m, 3);	// last deter bit * (-1)^(3+1)	// calc matrix determinanat	deter = temp[0] + temp[1] + temp[2] + temp[3];	// if there is no determinanat for this matrix	if( deter == 0.0f )		return false;	/* Calculating the inverse of a matrix */	/* 	 * To calculate the inverse of a matrix we have to split	 * the 4x4 into submatrices for each element (ignore row and column).	 * Then we take the determinant of the submatrix, and divide that by	 * the determinant of the whole matrix. Then muliply by (-1)^(r+c)	 * where r is the row and c is the column. NOTE: both of these should	 * start at 1, not 0.	 */	// calculate the inverse bit for the top row of the matrix (what	// we just calculated to find the determinant of the 4x4)	for(int i=0; i < 4; i++){		temp /= deter;		// For the first row, the row is always the same (1), and		// only the column changes.		// (i+1) to get it to one indexed instead of zero indexed,		// and the extra 1 because the row stays the same all the time.		// Also, (-1)^(even number) is always positive and		// (-1)^(odd number) is always negative.		if( (i+1)+1 % 2 )			temp = -temp;	}	int el = 0;	// now split the rest of the matrix into 3x3 submatrices	for(i=1; i < 4; i++){		for(int j=0; j < 4; j++){			el = i*4 + j;			split(m, j, i);			temp[el] = calcDeter(m, 3) / deter;						// negify is needed (yeah, that's probably not a real word)			if( (i+1)+(j+1) % 2 )				temp[el] = -temp[el];		}// for j < 4	}// 	// copy across	memcpy(m_mat, temp, sizeof(float)*16);	(*this).transpose();	// to finally get the inverse	// FIXME - not finished yet 3/01/06	return true;	// the inverse was successful}/** * Calculate the deteminant of a 3x3 or 2x2 matrix * \param m A pointer to the matrix values in an array * \param size The size of the matrix (can only be a square matrix) */float CMatrix::calcDeter(const float* m, int size) const{	// calculate the determinant of a 2x2	if( size == 2 ){		return (m[0] * m[1]) - (m[2] * m[3]);	}	// calculate the deteminant of a 3x3	else if( size == 3 ){		float d[3];	// hold the determinant at each stage		float a[4];	// the temp array to hold each 2x2 submatrix		// NOTE: in maths the first column of an array is column 1		// NOT 0. It is only 0 because arrays in C are zero-indexed,		// whereas matrix math appears to be one indexed.		// block top and column 1, so (deter * (-1)^(1+1))		a[0] = m[4];	a[1] = m[5];		a[2] = m[7];	a[3] = m[8];		d[0] = calcDeter(a,2);		// block top and column 2, so (deter * (-1)^(2+1))		a[0] = m[3];	// don't have to change		a[2] = m[6];	// el 1 or el 3		d[1] = -calcDeter(a,2);		// block top and column 3, so (deter * (-1)^(3+1))		a[1] = m[4];	// don't have to change		a[3] = m[7];	// el 0 or el 3		d[2] = calcDeter(a,2);		// calc the determinant of the 3x3 submatrix		return d[0] + d[1] + d[2];	}// if size == 3	// problem. Given an invalid size return 0 as the deteminant	return 0;}/** * Split 'this' 4x4 matrix into a smaller 3x3 submatrix depending on * which row and column to ignore. (used in the inverse function) * \param m A pointer to a float array to place the 3x3 submatrix (needs 9 elements) * \param colig The column to ignore * \param rowig The row to ignore */void CMatrix::split( float* m, int colig, int rowig) const{	int arr_count = 0;		// keep track of what element 					// in the array we're accessing	// split 'this' 4x4 into a 3x3	for(int i=0; i < 4; i++){		if( rowig == i )	// if we are ignoring this row			continue;	// skip back to the start of the loop		// traverse matrix column wise		for(int j=0; j < 4; j++){			if( colig == j )	// if we are ignoring this column				continue;	// skip back to start of the loop			// assign array value			m[ arr_count ] = (*this)(j,i);			arr_count++;	// only increment arr_count on element assignment		}// for j < 4	}// for i < 4}

// block top and column 1, so (deter * (-1)^(1+1))a[0] = m[4];	a[1] = m[5];a[2] = m[7];	a[3] = m[8];d[0] = calcDeter(a,2); // block top and column 2, so (deter * (-1)^(2+1))a[0] = m[3];	// don't have to changea[2] = m[6];	// el 1 or el 3d[1] = -calcDeter(a,2);

float& operator()(std::size_t i, std::size_t j){   assert(i < size && j < size);   return data[j*size + i];}const float& operator()(std::size_t i, std::size_t j) const{   assert(i < size && j < size);   return data[j*size + i];}

// block top and column 1, so (deter * (-1)^(1+1))a[0] = m[4];	a[1] = m[5];a[2] = m[7];	a[3] = m[8];d[0] = calcDeter(a,2); // block top and column 2, so (deter * (-1)^(2+1))a[0] = m[3];	// don't have to changea[2] = m[6];	// el 1 or el 3d[1] = -calcDeter(a,2);

split(m, 0, 0);temp[0] = calcDeter(m, 3);	// first deter bit * (-1)^((0+1)+1)split(m, 1, 0);temp[1] = -calcDeter(m, 3);	// second deter bit * (-1)^((1+1)+1)		split(m, 2, 0);temp[2] = calcDeter(m, 3);	// third deter bit * (-1)^((2+1)+1)split(m, 3, 0);temp[3] = -calcDeter(m, 3);	// last deter bit * (-1)^((3+1)+1)// calc matrix determinanatdeter = temp[0] + temp[1] + temp[2] + temp[3];