Aquí hay una versión de la respuesta de batty, pero esto calcula el correcto inverso. la versión de Batty calcula la transposición de la inversa.

// computes the inverse of a matrix m
double det = m(0, 0) * (m(1, 1) * m(2, 2) - m(2, 1) * m(1, 2)) -
             m(0, 1) * (m(1, 0) * m(2, 2) - m(1, 2) * m(2, 0)) +
             m(0, 2) * (m(1, 0) * m(2, 1) - m(1, 1) * m(2, 0));

double invdet = 1 / det;

Matrix33d minv; // inverse of matrix m
minv(0, 0) = (m(1, 1) * m(2, 2) - m(2, 1) * m(1, 2)) * invdet;
minv(0, 1) = (m(0, 2) * m(2, 1) - m(0, 1) * m(2, 2)) * invdet;
minv(0, 2) = (m(0, 1) * m(1, 2) - m(0, 2) * m(1, 1)) * invdet;
minv(1, 0) = (m(1, 2) * m(2, 0) - m(1, 0) * m(2, 2)) * invdet;
minv(1, 1) = (m(0, 0) * m(2, 2) - m(0, 2) * m(2, 0)) * invdet;
minv(1, 2) = (m(1, 0) * m(0, 2) - m(0, 0) * m(1, 2)) * invdet;
minv(2, 0) = (m(1, 0) * m(2, 1) - m(2, 0) * m(1, 1)) * invdet;
minv(2, 1) = (m(2, 0) * m(0, 1) - m(0, 0) * m(2, 1)) * invdet;
minv(2, 2) = (m(0, 0) * m(1, 1) - m(1, 0) * m(0, 1)) * invdet;

Este fragmento de código calcula el inverso transpuesto de la matriz A:

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;

Aunque la pregunta estipulaba matrices no singulares, es posible que aún desee verificar si el determinante es igual a cero (o muy cerca de cero) y marcarlo de alguna manera para estar seguro.

Con el debido respeto a nuestro desconocido (yahoo) póster, miro código así y simplemente muero un poco por dentro. La sopa de letras es increíblemente difícil de depurar. Un solo error tipográfico en cualquier lugar puede arruinar todo el día. Lamentablemente, este ejemplo en particular carecía de variables con guiones bajos. Es mucho más divertido cuando tenemos a_b-c_d * e_f-g_h. Especialmente cuando se usa una fuente donde _ y-tienen la misma longitud de píxel.

Retomando Suvesh Pratapa sobre su sugerencia, nota:

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

(A) Al tomar un menor de una matriz 3x3, tenemos 4 valores de interés. El índice X/Y más bajo es siempre 0 o 1. El índice X/Y más alto es siempre 1 o 2. Siempre! Por lo tanto:

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) El determinante es ahora: (¡Note el signo menos!)

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) Y el inverso es ahora:

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

Y redondearlo con un código de prueba de menor calidad:

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

Pensado de último momento:

También es posible que desee detectar determinantes muy grandes como redondeo ¡los errores afectarán su precisión!

Un archivo de encabezado bastante agradable (creo) que contiene macros para la mayoría de las operaciones de matriz 2x2, 3x3 y 4x4 ha estado disponible con la mayoría de los kits de herramientas de OpenGL. No como estándar, pero lo he visto en varios lugares.

Puede comprobarlo aquí. Al final de la misma encontrará tanto inversa de 2x2, 3x3 y 4x4.

Vvector.h

No trate de hacer esto usted mismo si usted es serio acerca de conseguir casos edge derecho. Así que mientras que muchos métodos ingenuos / simples son teóricamente exactos, pueden tener un comportamiento numérico desagradable para matrices casi singulares. En particular, puede obtener errores de cancelación/redondeo que hacen que obtenga resultados arbitrariamente malos.

Una forma "correcta" es la eliminación gaussiana con pivote de fila y columna para que siempre esté dividiendo por el mayor valor numérico restante. (Esto también es estable para matrices NxN.). Tenga en cuenta que el pivote de fila por sí solo no captura todos los casos malos.

Sin embargo IMO implementando este derecho y rápido no vale la pena su tiempo-utilice una biblioteca bien probada y hay un montón de encabezados solo unos.

Acabo de crear una clase QMatrix. Utiliza el vector > contenedor incorporado. QMatrix.h Utiliza el método Jordan-Gauss para calcular la inversa de una matriz cuadrada.

Puede usarlo de la siguiente manera:

#include "QMatrix.h"
#include <iostream>

int main(){
QMatrix<double> A(3,3,true);
QMatrix<double> Result = A.inverse()*A; //should give the idendity matrix

std::cout<<A.inverse()<<std::endl;
std::cout<<Result<<std::endl; // for checking
return 0;
}

La función inversa se implementa de la siguiente manera:

Dada una clase con los siguientes campos:

template<class T> class QMatrix{
public:
int rows, cols;
std::vector<std::vector<T> > A;

La función inversa ():

template<class T> 
QMatrix<T> QMatrix<T>:: inverse(){
Identity<T> Id(rows); //the Identity Matrix as a subclass of QMatrix.
QMatrix<T> Result = *this; // making a copy and transforming it to the Identity matrix
T epsilon = 0.000001;
for(int i=0;i<rows;++i){
    //check if Result(i,i)==0, if true, switch the row with another

    for(int j=i;j<rows;++j){
        if(std::abs(Result(j,j))<epsilon) { //uses Overloading()(int int) to extract element from Result Matrix
            Result.replace_rows(i,j+1); //switches rows i with j+1
        }
        else break;
    }
    // main part, making a triangular matrix
    Id(i)=Id(i)*(1.0/Result(i,i));
    Result(i)=Result(i)*(1.0/Result(i,i));  // Using overloading ()(int) to get a row form the matrix
    for(int j=i+1;j<rows;++j){
        T temp = Result(j,i);
        Result(j) = Result(j) - Result(i)*temp;
        Id(j) = Id(j) - Id(i)*temp; //doing the same operations to the identity matrix
        Result(j,i)=0; //not necessary, but looks nicer than 10^-15
    }
}

// solving a triangular matrix 
for(int i=rows-1;i>0;--i){
    for(int j=i-1;j>=0;--j){
        T temp = Result(j,i);
        Id(j) = Id(j) - temp*Id(i);
        Result(j)=Result(j)-temp*Result(i);
    }
}

return Id;
}

También recomendaría Ilmbase, que es parte de OpenEXR. Es un buen conjunto de plantillas 2,3,4-vector y rutinas de matriz.

# 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 += 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;
    } /* end do */

    while (ch == 'y');
    getch ();

    return 0;
}

#include <iostream>
using namespace std;

int main()
{
    double A11, A12, A13;
    double A21, A22, A23;
    double A31, A32, A33;

    double B11, B12, B13;
    double B21, B22, B23;
    double B31, B32, B33;

    cout << "Enter all number from left to right, from top to bottom, and press enter after every number: ";
    cin  >> A11;
    cin  >> A12;
    cin  >> A13;
    cin  >> A21;
    cin  >> A22;
    cin  >> A23;
    cin  >> A31;
    cin  >> A32;
    cin  >> A33;

    B11 = 1 / ((A22 * A33) - (A23 * A32));
    B12 = 1 / ((A13 * A32) - (A12 * A33));
    B13 = 1 / ((A12 * A23) - (A13 * A22));
    B21 = 1 / ((A23 * A31) - (A21 * A33));
    B22 = 1 / ((A11 * A33) - (A13 * A31));
    B23 = 1 / ((A13 * A21) - (A11 * A23));
    B31 = 1 / ((A21 * A32) - (A22 * A31));
    B32 = 1 / ((A12 * A31) - (A11 * A32));
    B33 = 1 / ((A11 * A22) - (A12 * A21));

    cout << B11 << "\t" << B12 << "\t" << B13 << endl;
    cout << B21 << "\t" << B22 << "\t" << B23 << endl;
    cout << B31 << "\t" << B32 << "\t" << B33 << endl;

    return 0;
}

//Title: Matrix Header File
//Writer: Say OL
//This is a beginner code not an expert one
//No responsibilty for any errors
//Use for your own risk
using namespace std;
int row,col,Row,Col;
double Coefficient;
//Input Matrix
void Input(double Matrix[9][9],int Row,int Col)
{
    for(row=1;row<=Row;row++)
        for(col=1;col<=Col;col++)
        {
            cout<<"e["<<row<<"]["<<col<<"]=";
            cin>>Matrix[row][col];
        }
}
//Output Matrix
void Output(double Matrix[9][9],int Row,int Col)
{
    for(row=1;row<=Row;row++)
    {
        for(col=1;col<=Col;col++)
            cout<<Matrix[row][col]<<"\t";
        cout<<endl;
    }
}
//Copy Pointer to Matrix
void CopyPointer(double (*Pointer)[9],double Matrix[9][9],int Row,int Col)
{
    for(row=1;row<=Row;row++)
        for(col=1;col<=Col;col++)
            Matrix[row][col]=Pointer[row][col];
}
//Copy Matrix to Matrix
void CopyMatrix(double MatrixInput[9][9],double MatrixTarget[9][9],int Row,int Col)
{
    for(row=1;row<=Row;row++)
        for(col=1;col<=Col;col++)
            MatrixTarget[row][col]=MatrixInput[row][col];
}
//Transpose of Matrix
double MatrixTran[9][9];
double (*(Transpose)(double MatrixInput[9][9],int Row,int Col))[9]
{
    for(row=1;row<=Row;row++)
        for(col=1;col<=Col;col++)
            MatrixTran[col][row]=MatrixInput[row][col];
    return MatrixTran;
}
//Matrix Addition
double MatrixAdd[9][9];
double (*(Addition)(double MatrixA[9][9],double MatrixB[9][9],int Row,int Col))[9]
{
    for(row=1;row<=Row;row++)
        for(col=1;col<=Col;col++)
            MatrixAdd[row][col]=MatrixA[row][col]+MatrixB[row][col];
    return MatrixAdd;
}
//Matrix Subtraction
double MatrixSub[9][9];
double (*(Subtraction)(double MatrixA[9][9],double MatrixB[9][9],int Row,int Col))[9]
{
    for(row=1;row<=Row;row++)
        for(col=1;col<=Col;col++)
            MatrixSub[row][col]=MatrixA[row][col]-MatrixB[row][col];
    return MatrixSub;
}
//Matrix Multiplication
int mRow,nCol,pCol,kcol;
double MatrixMult[9][9];
double (*(Multiplication)(double MatrixA[9][9],double MatrixB[9][9],int mRow,int nCol,int pCol))[9]
{
    for(row=1;row<=mRow;row++)
        for(col=1;col<=pCol;col++)
        {
            MatrixMult[row][col]=0.0;
            for(kcol=1;kcol<=nCol;kcol++)
                MatrixMult[row][col]+=MatrixA[row][kcol]*MatrixB[kcol][col];
        }
    return MatrixMult;
}
//Interchange Two Rows
double RowTemp[9][9];
double MatrixInter[9][9];
double (*(InterchangeRow)(double MatrixInput[9][9],int Row,int Col,int iRow,int jRow))[9]
{
    CopyMatrix(MatrixInput,MatrixInter,Row,Col);
    for(col=1;col<=Col;col++)
    {
        RowTemp[iRow][col]=MatrixInter[iRow][col];
        MatrixInter[iRow][col]=MatrixInter[jRow][col];
        MatrixInter[jRow][col]=RowTemp[iRow][col];
    }
    return MatrixInter;
}
//Pivote Downward
double MatrixDown[9][9];
double (*(PivoteDown)(double MatrixInput[9][9],int Row,int Col,int tRow,int tCol))[9]
{
    CopyMatrix(MatrixInput,MatrixDown,Row,Col);
    Coefficient=MatrixDown[tRow][tCol];
    if(Coefficient!=1.0)
        for(col=1;col<=Col;col++)
            MatrixDown[tRow][col]/=Coefficient;
    if(tRow<Row)
        for(row=tRow+1;row<=Row;row++)
        {
            Coefficient=MatrixDown[row][tCol];
            for(col=1;col<=Col;col++)
                MatrixDown[row][col]-=Coefficient*MatrixDown[tRow][col];
        }
return MatrixDown;
}
//Pivote Upward
double MatrixUp[9][9];
double (*(PivoteUp)(double MatrixInput[9][9],int Row,int Col,int tRow,int tCol))[9]
{
    CopyMatrix(MatrixInput,MatrixUp,Row,Col);
    Coefficient=MatrixUp[tRow][tCol];
    if(Coefficient!=1.0)
        for(col=1;col<=Col;col++)
            MatrixUp[tRow][col]/=Coefficient;
    if(tRow>1)
        for(row=tRow-1;row>=1;row--)
        {
            Coefficient=MatrixUp[row][tCol];
            for(col=1;col<=Col;col++)
                MatrixUp[row][col]-=Coefficient*MatrixUp[tRow][col];
        }
    return MatrixUp;
}
//Pivote in Determinant
double MatrixPiv[9][9];
double (*(Pivote)(double MatrixInput[9][9],int Dim,int pTarget))[9]
{
    CopyMatrix(MatrixInput,MatrixPiv,Dim,Dim);
    for(row=pTarget+1;row<=Dim;row++)
    {
        Coefficient=MatrixPiv[row][pTarget]/MatrixPiv[pTarget][pTarget];
        for(col=1;col<=Dim;col++)
        {
            MatrixPiv[row][col]-=Coefficient*MatrixPiv[pTarget][col];
        }
    }
    return MatrixPiv;
}
//Determinant of Square Matrix
int dCounter,dRow;
double Det;
double MatrixDet[9][9];
double Determinant(double MatrixInput[9][9],int Dim)
{
    CopyMatrix(MatrixInput,MatrixDet,Dim,Dim);
    Det=1.0;
    if(Dim>1)
    {
        for(dRow=1;dRow<Dim;dRow++)
        {
            dCounter=dRow;
            while((MatrixDet[dRow][dRow]==0.0)&(dCounter<=Dim))
            {
                dCounter++;
                Det*=-1.0;
                CopyPointer(InterchangeRow(MatrixDet,Dim,Dim,dRow,dCounter),MatrixDet,Dim,Dim);
            }
            if(MatrixDet[dRow][dRow]==0)
            {
                Det=0.0;
                break;
            }
            else
            {
                Det*=MatrixDet[dRow][dRow];
                CopyPointer(Pivote(MatrixDet,Dim,dRow),MatrixDet,Dim,Dim);
            }
        }
        Det*=MatrixDet[Dim][Dim];
    }
    else Det=MatrixDet[1][1];
    return Det;
}
//Matrix Identity
double MatrixIdent[9][9];
double (*(Identity)(int Dim))[9]
{
    for(row=1;row<=Dim;row++)
        for(col=1;col<=Dim;col++)
            if(row==col)
                MatrixIdent[row][col]=1.0;
            else
                MatrixIdent[row][col]=0.0;
    return MatrixIdent;
}
//Join Matrix to be Augmented Matrix
double MatrixJoin[9][9];
double (*(JoinMatrix)(double MatrixA[9][9],double MatrixB[9][9],int Row,int ColA,int ColB))[9]
{
    Col=ColA+ColB;
    for(row=1;row<=Row;row++)
        for(col=1;col<=Col;col++)
            if(col<=ColA)
                MatrixJoin[row][col]=MatrixA[row][col];
            else
                MatrixJoin[row][col]=MatrixB[row][col-ColA];
    return MatrixJoin;
}
//Inverse of Matrix
double (*Pointer)[9];
double IdentMatrix[9][9];
int Counter;
double MatrixAug[9][9];
double MatrixInv[9][9];
double (*(Inverse)(double MatrixInput[9][9],int Dim))[9]
{
    Row=Dim;
    Col=Dim+Dim;
    Pointer=Identity(Dim);
    CopyPointer(Pointer,IdentMatrix,Dim,Dim);
    Pointer=JoinMatrix(MatrixInput,IdentMatrix,Dim,Dim,Dim);
    CopyPointer(Pointer,MatrixAug,Row,Col);
    for(Counter=1;Counter<=Dim;Counter++)   
    {
        Pointer=PivoteDown(MatrixAug,Row,Col,Counter,Counter);
        CopyPointer(Pointer,MatrixAug,Row,Col);
    }
    for(Counter=Dim;Counter>1;Counter--)
    {
        Pointer=PivoteUp(MatrixAug,Row,Col,Counter,Counter);
        CopyPointer(Pointer,MatrixAug,Row,Col);
    }
    for(row=1;row<=Dim;row++)
        for(col=1;col<=Dim;col++)
            MatrixInv[row][col]=MatrixAug[row][col+Dim];
    return MatrixInv;
}
//Gauss-Jordan Elemination
double MatrixGJ[9][9];
double VectorGJ[9][9];
double (*(GaussJordan)(double MatrixInput[9][9],double VectorInput[9][9],int Dim))[9]
{
    Row=Dim;
    Col=Dim+1;
    Pointer=JoinMatrix(MatrixInput,VectorInput,Dim,Dim,1);
    CopyPointer(Pointer,MatrixGJ,Row,Col);
    for(Counter=1;Counter<=Dim;Counter++)   
    {
        Pointer=PivoteDown(MatrixGJ,Row,Col,Counter,Counter);
        CopyPointer(Pointer,MatrixGJ,Row,Col);
    }
    for(Counter=Dim;Counter>1;Counter--)
    {
        Pointer=PivoteUp(MatrixGJ,Row,Col,Counter,Counter);
        CopyPointer(Pointer,MatrixGJ,Row,Col);
    }
    for(row=1;row<=Dim;row++)
        for(col=1;col<=1;col++)
            VectorGJ[row][col]=MatrixGJ[row][col+Dim];
    return VectorGJ;
}
//Generalized Gauss-Jordan Elemination
double MatrixGGJ[9][9];
double VectorGGJ[9][9];
double (*(GeneralizedGaussJordan)(double MatrixInput[9][9],double VectorInput[9][9],int Dim,int vCol))[9]
{
    Row=Dim;
    Col=Dim+vCol;
    Pointer=JoinMatrix(MatrixInput,VectorInput,Dim,Dim,vCol);
    CopyPointer(Pointer,MatrixGGJ,Row,Col);
    for(Counter=1;Counter<=Dim;Counter++)   
    {
        Pointer=PivoteDown(MatrixGGJ,Row,Col,Counter,Counter);
        CopyPointer(Pointer,MatrixGGJ,Row,Col);
    }
    for(Counter=Dim;Counter>1;Counter--)
    {
        Pointer=PivoteUp(MatrixGGJ,Row,Col,Counter,Counter);
        CopyPointer(Pointer,MatrixGGJ,Row,Col);
    }
    for(row=1;row<=Row;row++)
        for(col=1;col<=vCol;col++)
            VectorGGJ[row][col]=MatrixGGJ[row][col+Dim];
    return VectorGGJ;
}
//Matrix Sparse, Three Diagonal Non-Zero Elements
double MatrixSpa[9][9];
double (*(Sparse)(int Dimension,double FirstElement,double SecondElement,double ThirdElement))[9]
{
    MatrixSpa[1][1]=SecondElement;
    MatrixSpa[1][2]=ThirdElement;
    MatrixSpa[Dimension][Dimension-1]=FirstElement;
    MatrixSpa[Dimension][Dimension]=SecondElement;
    for(int Counter=2;Counter<Dimension;Counter++)
    {
        MatrixSpa[Counter][Counter-1]=FirstElement;
        MatrixSpa[Counter][Counter]=SecondElement;
        MatrixSpa[Counter][Counter+1]=ThirdElement;
    }
    return MatrixSpa;
}

Copie y guarde el código anterior como Matriz.h a continuación, intente el siguiente código:

#include<iostream>
#include<conio.h>
#include"Matrix.h"
int Dim;
double Matrix[9][9];
int main()
{
    cout<<"Enter your matrix dimension: ";
    cin>>Dim;
    Input(Matrix,Dim,Dim);
    cout<<"Your matrix:"<<endl;
    Output(Matrix,Dim,Dim);
    cout<<"The inverse:"<<endl;
    Output(Inverse(Matrix,Dim),Dim,Dim);
    getch();
}

//Function for inverse of the input square matrix 'J' of dimension 'dim':

vector<vector<double > > inverseVec33(vector<vector<double > > J, int dim)
{
//Matrix of Minors
 vector<vector<double > > invJ(dim,vector<double > (dim));
for(int i=0; i<dim; i++)
{
    for(int j=0; j<dim; j++)
    {
        invJ[i][j] = (J[(i+1)%dim][(j+1)%dim]*J[(i+2)%dim][(j+2)%dim] -
                      J[(i+2)%dim][(j+1)%dim]*J[(i+1)%dim][(j+2)%dim]);
    }
}

//determinant of the matrix:
double detJ = 0.0;
for(int j=0; j<dim; j++)
{ detJ += J[0][j]*invJ[0][j];}

//Inverse of the given matrix.
 vector<vector<double > > invJT(dim,vector<double > (dim));
 for(int i=0; i<dim; i++)
{
    for(int j=0; j<dim; j++)
    {
        invJT[i][j] = invJ[j][i]/detJ;
    }
}

return invJT;
}

void main()
{
    //given matrix:
vector<vector<double > > Jac(3,vector<double > (3));
Jac[0][0] = 1; Jac[0][1] = 2;  Jac[0][2] = 6;
Jac[1][0] = -3; Jac[1][1] = 4;  Jac[1][2] = 3;
Jac[2][0] = 5; Jac[2][1] = 1;  Jac[2][2] = -4;`

//Inverse of the matrix Jac:
vector<vector<double > > JacI(3,vector<double > (3));
    //call function and store inverse of J as JacI:
JacI = inverseVec33(Jac,3);
}

Seguí adelante y lo escribí en python ya que creo que es mucho más legible que en c++ para un problema como este. El orden de la función está en orden de operaciones para resolver esto a mano a través de este video. Solo importa esto y llama a "print_invert" en tu matriz.

def print_invert (matrix):
  i_matrix = invert_matrix (matrix)
  for line in i_matrix:
    print (line)
  return

def invert_matrix (matrix):
  determinant = str (determinant_of_3x3 (matrix))
  cofactor = make_cofactor (matrix)
  trans_matrix = transpose_cofactor (cofactor)

  trans_matrix[:] = [[str (element) +'/'+ determinant for element in row] for row in trans_matrix]

  return trans_matrix

def determinant_of_3x3 (matrix):
  multiplication = 1
  neg_multiplication = 1
  total = 0
  for start_column in range (3):
    for row in range (3):
      multiplication *= matrix[row][(start_column+row)%3]
      neg_multiplication *= matrix[row][(start_column-row)%3]
    total += multiplication - neg_multiplication
    multiplication = neg_multiplication = 1
  if total == 0:
    total = 1
  return total

def make_cofactor (matrix):
  cofactor = [[0,0,0],[0,0,0],[0,0,0]]
  matrix_2x2 = [[0,0],[0,0]]
  # For each element in matrix...
  for row in range (3):
    for column in range (3):

      # ...make the 2x2 matrix in this inner loop
      matrix_2x2 = make_2x2_from_spot_in_3x3 (row, column, matrix)
      cofactor[row][column] = determinant_of_2x2 (matrix_2x2)

  return flip_signs (cofactor)

def make_2x2_from_spot_in_3x3 (row, column, matrix):
  c_count = 0
  r_count = 0
  matrix_2x2 = [[0,0],[0,0]]
  # ...make the 2x2 matrix in this inner loop
  for inner_row in range (3):
    for inner_column in range (3):
      if row is not inner_row and inner_column is not column:
        matrix_2x2[r_count % 2][c_count % 2] = matrix[inner_row][inner_column]
        c_count += 1
    if row is not inner_row:
      r_count += 1
  return matrix_2x2

def determinant_of_2x2 (matrix):
  total = matrix[0][0] * matrix [1][1]
  return total - (matrix [1][0] * matrix [0][1])

def flip_signs (cofactor):
  sign_pos = True 
  # For each element in matrix...
  for row in range (3):
    for column in range (3):
      if sign_pos:
        sign_pos = False
      else:
        cofactor[row][column] *= -1
        sign_pos = True
  return cofactor

def transpose_cofactor (cofactor):
  new_cofactor = [[0,0,0],[0,0,0],[0,0,0]]
  for row in range (3):
    for column in range (3):
      new_cofactor[column][row] = cofactor[row][column]
  return new_cofactor