Help: 4x4 Matrix Inverse Implementation

[zb0S216C]

I'm implementing a 4x4 matrix class and all is going well until the inverse function turned up. I'm trying to implement the inverse function, but I can't seem to get my head around it.

I've tried the internet, but found nothing useful. Also, I've looked into source code of other programs/libraries that implement a matrix class, but the code is unreadable.

Any general idea how I can implement this damn 4x4 inverse function? I know the process of inversion, but putting that process into code is proving quite a challenge.

In addition, I do have some code, but it's unmanageable and inefficient at the moment. If you wan't to see it, just ask.

Additional question(s): What applications does the inverse matrix have in 3-D?

Thanks.

_Wazzak

[ne555]

>I do have some code, but it's unmanageable and inefficient at the moment.
> If you want to see it, just ask.
it would be nice to have something to work with.
¿inefficient? It's a 4x4 matrix, ¿how much could you screw up?

Using Gauss-Jordan and assuming that you do not need to change pivoting

matrix inverse(matrix M){
   matrix I = identity(M.size());
   for(int K=0; K<rows; ++K){
      //one in the pivot
      double factor = M[K][K];
      M[K] /= factor;
      I[K] /= factor;
      //zeroing the column
      for(int L=0; L<rows; ++L){
         if( K==L ) continue;
         double coefficient = M[L][K];
         M[L] -= coefficient*M[K];
         I[L] -= coefficient*I[K];
      }
   }

   return I;
}

[mik2718]

I am using a routine lifted from numerical recipes (2nd ed)

This edition is free to view on line here
http://apps.nrbook.com/c/index.html

This routine uses full pivoting and I have found it reliable. It will probably be inaccurate for near singular matrices. You need to use SVD and then least squares for that (there is also a routine for this in the book)

If you are doing geometric work then you can avoid singular matrices by considering the geometry of the situation.

I don't think routines without at least partial pivoting are useful. So its not that simple a problem.

Here's my code, no I don't understand every detail, it just works.

NB: You ask about inverses in 3D, Clearly if a 3x3 matrix represents a transformation of an object then the inverse will take it back to its original position. BUT when you are dealing with rotations (often the case) these are unit matrices, so the inverse is just the transpose (MUCH quicker to calculate). So unless you are doing unusual tranformations, code for the general inverse is not required.

class Mat4
{
  double **ptr;  // initialised to 4x4 array
public :
  void inverse(); // in place inverse
};

void Mat4::inverse()
{ // full pivoting routine nicked from numerical recipes
  int indxc[4],indxr[4],ipiv[4];
  int i,icol=0,irow=0,j,k,l,ll;
  double big,dum,pivinv,temp;
  double* temp2;
  for(j=0;j<4;++j)
    ipiv[j] = 0;
  for(i=0;i<4;++i)
  {
    big = 0.0;
    for(j=0;j<4;++j)
    {
      if(ipiv[j] != 1)
      {
        for(k=0;k<4;++k)
	{
	  if(ipiv[k] == 0)
	  {
	    if(fabs(ptr[j][k]) >= big)
	    {
	      big = fabs(ptr[j][k]);
	      irow = j;
	      icol = k;
	    }
	  }
	  else
	  {
	    if(ipiv[k] > 1)
	    {
              throw Inv4DSingular();
	    }
	  }
	}
      }
    }
    ++ipiv[icol];
    if(irow != icol)
    {
      temp2 = ptr[irow];
      ptr[irow] = ptr[icol];
      ptr[icol] = temp2;
    }
    indxr[i] = irow;
    indxc[i] = icol;
    if(ptr[icol][icol] == 0.0)
    {
      throw Inv4DSingular();
    }
    pivinv = 1.0/ptr[icol][icol];
    ptr[icol][icol] = 1.0;
    for(l=0;l<4;++l)
      ptr[icol][l] *= pivinv;
    for(ll=0;ll<4;++ll)
    {
      if(ll != icol)
      {
        dum = ptr[ll][icol];
	ptr[ll][icol] = 0.0;
	for(l=0;l<4;++l)
	  ptr[ll][l] -= ptr[icol][l]*dum;
      }
    }
  }
  for(l=3;l>=0;--l)
  {
    if(indxr[l] != indxc[l])
    {
      for(k=0;k<4;++k)
      {
        temp = ptr[k][indxr[l]];
	ptr[k][indxr[l]] = ptr[k][indxc[l]];
	ptr[k][indxc[l]] = temp;
      }
    }
  }
}

[zb0S216C]

@ne555: Thanks for the response. Using Gauss-Jordon in a way that's efficient is proving quite the challenge.

@mik2718: Thanks for the response. The implementation you provided is a little too bloated for my needs. That exception, and all of those loops and conditionals are off-putting. Thanks, though :)

Update:

I've decided to push aside my version of the 4x4 inversion function in favour of GLM's matrix inversion function. It's relatively neat and branch-less; I've converted it to fit my matrix class.

Thanks for responding :)

_Wazzak