Simple 3x3 matrix inverse code (C++)

Question

What's the easiest way to compute a 3x3 matrix inverse?

I'm just looking for a short code snippet that'll do the trick for non-singular matrices, possibly using Cramer's rule. It doesn't need to be highly optimized. I'd prefer simplicity over speed. I'd rather not link in additional libraries. Primarily I was hoping to have this on Stack Overflow so that I wouldn't have to hunt around for it or rewrite from scratch again next time.

Answer 1

Look here what google helped me find ! :)

Answer 2

You give no details and your question is very generic, so my answer is to use BLAS.

Answer 3

This seems to work:

double determinant =    +A(0,0)*(A(1,1)*A(2,2)-A(2,1)*A(1,2))
                        -A(0,1)*(A(1,0)*A(2,2)-A(1,2)*A(2,0))
                        +A(0,2)*(A(1,0)*A(2,1)-A(1,1)*A(2,0));
double invdet = 1/determinant;
result(0,0) =  (A(1,1)*A(2,2)-A(2,1)*A(1,2))*invdet;
result(1,0) = -(A(0,1)*A(2,2)-A(0,2)*A(2,1))*invdet;
result(2,0) =  (A(0,1)*A(1,2)-A(0,2)*A(1,1))*invdet;
result(0,1) = -(A(1,0)*A(2,2)-A(1,2)*A(2,0))*invdet;
result(1,1) =  (A(0,0)*A(2,2)-A(0,2)*A(2,0))*invdet;
result(2,1) = -(A(0,0)*A(1,2)-A(1,0)*A(0,2))*invdet;
result(0,2) =  (A(1,0)*A(2,1)-A(2,0)*A(1,1))*invdet;
result(1,2) = -(A(0,0)*A(2,1)-A(2,0)*A(0,1))*invdet;
result(2,2) =  (A(0,0)*A(1,1)-A(1,0)*A(0,1))*invdet;

Though the question stipulated non-singular matrices, you might still want to check if determinant equals zero (or very near zero) and flag it in some way to be safe.

Answer 4

With all due respect to our unknown (yahoo) poster, I look at code like that and just die a little inside. Alphabet soup is just so insanely difficult to debug. A single typo anywhere in there can really ruin your whole day. Sadly, this particular example lacked variables with underscores. It's so much more fun when we have a_b-c_d*e_f-g_h. Especially when using a font where _ and - have the same pixel length.

Taking up Suvesh Pratapa on his suggestion, I note:

Given 3x3 matrix:
       y0x0  y0x1  y0x2
       y1x0  y1x1  y1x2
       y2x0  y2x1  y2x2
Declared as double matrix [/*Y=*/3] [/*X=*/3];

(A) When taking a minor of a 3x3 array, we have 4 values of interest. The lower X/Y index is always 0 or 1. The higher X/Y index is always 1 or 2. Always! Therefore:

double determinantOfMinor( int          theRowHeightY,
                           int          theColumnWidthX,
                           const double theMatrix [/*Y=*/3] [/*X=*/3] )
{
  int x1 = theColumnWidthX == 0 ? 1 : 0;  /* always either 0 or 1 */
  int x2 = theColumnWidthX == 2 ? 1 : 2;  /* always either 1 or 2 */
  int y1 = theRowHeightY   == 0 ? 1 : 0;  /* always either 0 or 1 */
  int y2 = theRowHeightY   == 2 ? 1 : 2;  /* always either 1 or 2 */

  return ( theMatrix [y1] [x1]  *  theMatrix [y2] [x2] )
      -  ( theMatrix [y1] [x2]  *  theMatrix [y2] [x1] );
}

(B) Determinant is now: (Note the minus sign!)

double determinant( const double theMatrix [/*Y=*/3] [/*X=*/3] )
{
  return ( theMatrix [0] [0]  *  determinantOfMinor( 0, 0, theMatrix ) )
      -  ( theMatrix [0] [1]  *  determinantOfMinor( 0, 1, theMatrix ) )
      +  ( theMatrix [0] [2]  *  determinantOfMinor( 0, 2, theMatrix ) );
}

(C) And the inverse is now:

bool inverse( const double theMatrix [/*Y=*/3] [/*X=*/3],
                    double theOutput [/*Y=*/3] [/*X=*/3] )
{
  double det = determinant( theMatrix );

    /* Arbitrary for now.  This should be something nicer... */
  if ( ABS(det) < 1e-2 )
  {
    memset( theOutput, 0, sizeof theOutput );
    return false;
  }

  double oneOverDeterminant = 1.0 / det;

  for (   int y = 0;  y < 3;  y ++ )
    for ( int x = 0;  x < 3;  x ++   )
    {
        /* Rule is inverse = 1/det * minor of the TRANSPOSE matrix.  *
         * Note (y,x) becomes (x,y) INTENTIONALLY here!              */
      theOutput [y] [x]
        = determinantOfMinor( x, y, theMatrix ) * oneOverDeterminant;

        /* (y0,x1)  (y1,x0)  (y1,x2)  and (y2,x1)  all need to be negated. */
      if( 1 == ((x + y) % 2) )
        theOutput [y] [x] = - theOutput [y] [x];
    }

  return true;
}

And round it out with a little lower-quality testing code:

void printMatrix( const double theMatrix [/*Y=*/3] [/*X=*/3] )
{
  for ( int y = 0;  y < 3;  y ++ )
  {
    cout << "[  ";
    for ( int x = 0;  x < 3;  x ++   )
      cout << theMatrix [y] [x] << "  ";
    cout << "]" << endl;
  }
  cout << endl;
}

void matrixMultiply(  const double theMatrixA [/*Y=*/3] [/*X=*/3],
                      const double theMatrixB [/*Y=*/3] [/*X=*/3],
                            double theOutput  [/*Y=*/3] [/*X=*/3]  )
{
  for (   int y = 0;  y < 3;  y ++ )
    for ( int x = 0;  x < 3;  x ++   )
    {
      theOutput [y] [x] = 0;
      for ( int i = 0;  i < 3;  i ++ )
        theOutput [y] [x] +=  theMatrixA [y] [i] * theMatrixB [i] [x];
    }
}

int
main(int argc, char **argv)
{
  if ( argc > 1 )
    SRANDOM( atoi( argv[1] ) );

  double m[3][3] = { { RANDOM_D(0,1e3), RANDOM_D(0,1e3), RANDOM_D(0,1e3) },
                     { RANDOM_D(0,1e3), RANDOM_D(0,1e3), RANDOM_D(0,1e3) },
                     { RANDOM_D(0,1e3), RANDOM_D(0,1e3), RANDOM_D(0,1e3) } };
  double o[3][3], mm[3][3];

  if ( argc <= 2 )
    cout << fixed << setprecision(3);

  printMatrix(m);
  cout << endl << endl;

  SHOW( determinant(m) );
  cout << endl << endl;

  BOUT( inverse(m, o) );
  printMatrix(m);
  printMatrix(o);
  cout << endl << endl;

  matrixMultiply (m, o, mm );
  printMatrix(m);
  printMatrix(o);
  printMatrix(mm);  
  cout << endl << endl;
}

---

Afterthought:

You may also want to detect very large determinants as round-off errors will affect your accuracy!

Answer 5

A rather nice (I think) header file containing macros for most 2x2, 3x3 and 4x4 matrix operations has been available with most OpenGL toolkits. Not as standard but I've seen it at various places.

You can check it out here. At the end of it you will find both inverse of 2x2, 3x3 and 4x4.

vvector.h

Answer 6

I would also recommend Ilmbase, which is part of OpenEXR. It's a good set of templated 2,3,4-vector and matrix routines.

Answer 7

# include <conio.h>
# include<iostream.h>
const int size=9;
int main()
{
char ch;
do
{
clrscr();
int i,j,x,y,z,det,a[size],b[size];
cout<<"           **** MATRIX OF 3x3 ORDER ****"<<endl<<endl<<endl;
for (i=0;i<=size;i++)
a[i]=0;
for (i=0;i<size;i++)
{
cout<<"Enter "<<i+1<<" element of matrix=";
cin>>a[i]; cout<<endl<<endl;
}
clrscr();
cout<<"your entered matrix is "<<endl<<endl;
for (i=0;i<size;i=i+3)
cout<<a[i]<<"  "<<a[i+1]<<"  "<<a[i+2]<<endl<<endl;
cout<<"Transpose of given matrix is"<<endl<<endl;
for (i=0;i<3;i++)
cout<<a[i]<<"  "<<a[i+3]<<"  "<<a[i+6]<<endl<<endl;
cout<<"Determinent of given matrix is=";
x=a[0]*(a[4]*a[8]-a[5]*a[7]);
y=a[1]*(a[3]*a[8]-a[5]*a[6]);
z=a[2]*(a[3]*a[7]-a[4]*a[6]);
det=x-y+z;
cout<<det<<endl<<endl<<endl<<endl;
if (det==0)
{
cout<<"As Determinent=0 so it is singular matrix and its inverse cannot exist"<<endl<<endl;
goto quit;
}
b[0]=a[4]*a[8]-a[5]*a[7];
b[1]=a[5]*a[6]-a[3]*a[8];
b[2]=a[3]*a[7]-a[4]*a[6];
b[3]=a[2]*a[7]-a[1]*a[8];
b[4]=a[0]*a[8]-a[2]*a[6];
b[5]=a[1]*a[6]-a[0]*a[7];
b[6]=a[1]*a[5]-a[2]*a[4];
b[7]=a[2]*a[3]-a[0]*a[5];
b[8]=a[0]*a[4]-a[1]*a[3];
cout<<"Adjoint of given matrix is"<<endl<<endl;
for (i=0;i<3;i++)
{
cout<<b[i]<<"  "<<b[i+3]<<"  "<<b[i+6]<<endl<<endl;
}
cout<<endl<<endl;
cout<<"Inverse of given matrix is "<<endl<<endl<<endl;
for (i=0;i<3;i++)
{
cout<<b[i]<<"/"<<det<<"  "<<b[i+3]<<"/"<<det<<"  "<<b[i+6]<<"/"<<det<<endl<<endl;
}
quit:
cout<<endl<<endl;
cout<<"Do You want to continue this again press (y/yes,n/no)";
cin>>ch; cout<<endl<<endl;
}
while (ch=='y');
getch ();
return 0;
}