简单的 3x3 矩阵逆码 (C++)

Question

计算 3x3 矩阵逆的最简单方法是什么？

我只是在寻找一个简短的代码片段，可以使用克莱默规则来解决非奇异矩阵的问题。 它不需要高度优化。 我更喜欢简单而不是速度。 我不想链接其他库。

Answer 1

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

Answer 2

我还推荐 Ilmbase，它是 OpenEXR 的一部分。 这是一组很好的模板化 2,3,4 向量和矩阵例程。

Answer 3

大多数 OpenGL 工具包都提供了一个相当不错的（我认为）头文件，其中包含大多数 2x2、3x3 和 4x4 矩阵运算的宏。 不是标准的，但我在不同的地方见过它。

你可以在这里查看。 最后你会发现 2x2、3x3 和 4x4 的逆。

vvector.h

Answer 4

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

Answer 5

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

复制上述代码并将其保存为 Matrix.h，然后尝试以下代码：

#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();
}

Answer 6

这是巴蒂答案的一个版本，但这计算出正确的逆。 batty 的版本计算逆矩阵的转置。

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

Answer 7

恕我直言，我对我们未知的（雅虎）发布者表示敬意，我看着这样的代码，心里有点死了。 字母汤实在是太难调试了。 任何地方的一个拼写错误都会毁掉你的一整天。 遗憾的是，这个特定的示例缺少带下划线的变量。 当我们有 a_b-c_d*e_f-g_h 时，会更有趣。 特别是当使用 _ 和 - 具有相同像素长度的字体时。

根据 Suvesh Pratapa 的建议，我注意到：

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

(A) 当取 3x3 数组的次要值时，我们有 4 个感兴趣的值。 较低的 X/Y 索引始终为 0 或 1。较高的 X/Y 索引始终为 1 或 2。始终！ 因此：

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) 行列式现在是：（注意减号！）

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) 倒数现在是：

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

用一些质量较低的测试代码将其四舍五入：

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

---

事后思考：

您可能还想检测非常大的行列式四舍五入误差会影响您的准确性！

Answer 8

你为什么不尝试自己编码呢？ 把它当作一个挑战。 :)

对于 3×3 矩阵

_{（来源：wolfram.com）}

矩阵的逆矩阵是

_{（来源：wolfram.com）}

我假设您知道矩阵 |A| 的行列式是什么 是。

图像 (c) Wolfram|Alpha 和
mathworld.wolfram (06-11-09,
22.06)

Answer 9

这段代码计算矩阵 A 的转置逆矩阵 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;

虽然问题规定了非奇异矩阵，但您可能仍然想检查行列式是否等于零（或非常接近零）并将其标记为某些确保安全的方法。

Answer 10

如果您真的想正确处理边缘情况，请不要尝试自己这样做。 因此，虽然许多朴素/简单的方法在理论上是精确的，但对于近乎奇异的矩阵，它们可能会产生令人讨厌的数值行为。 特别是，您可能会遇到取消/舍入错误，导致您得到任意糟糕的结果。

一种“正确”的方法是使用行和列旋转的高斯消除法，以便您始终除以最大的剩余数值。 （这对于 NxN 矩阵也是稳定的。）。 请注意，仅行旋转并不能捕获所有不良情况。

然而，在我看来，正确而快速地实现这一点并不值得你花时间——使用一个经过良好测试的库，并且有一堆只有标头的库。

Answer 11

我刚刚创建了一个 QMatrix 类。 它使用内置的向量> 容器。 QMatrix.h
它使用 Jordan-Gauss 方法来计算方阵的逆。

您可以按如下方式使用它：

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

inverse 函数的实现如下：

给定一个具有以下字段的类：

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

inverse() 函数：

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

Answer 12

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