C语言中的方阵求逆

Question

我写了一个n*n方阵的求逆函数。

void inverseMatrix(int n, float **matrix)
{
    float ratio,a;
    int i, j, k;

   for(i = 0; i < n; i++)
   {
      for(j = n; j < 2*n; j++)
      {
         if(i==(j-n))
            matrix[i][j] = 1.0;
         else
            matrix[i][j] = 0.0;
      }
  }

  for(i = 0; i < n; i++)
  {
      for(j = 0; j < n; j++)
      {
          if(i!=j)
          {
              ratio = matrix[j][i]/matrix[i][i];
              for(k = 0; k < 2*n; k++)
              {
                  matrix[j][k] -= ratio * matrix[i][k];
              }
          }
      }
  }

  for(i = 0; i < n; i++)
  {
      a = matrix[i][i];
      for(j = 0; j < 2*n; j++)
      {
          matrix[i][j] /= a;
      }
  }

//return matrix;
}

这在几乎所有情况下都工作得很好，但在某些情况下会失败，如下图所示：

1 1 1 0       1 1 1 0
1 1 2 0       1 1 2 0  
1 2 0 1       1 2 1 0
1 2 0 2       1 2 0 2

我忽略的情况可能是什么？

谢谢!

Answer 1

参见http://www.sourcecodesworld.com/source/show.asp?ScriptID=1086。

使用高斯乔丹算法

#include<stdio.h>
#include<stdlib.h>

int main()
{
    float **A,**I,temp;
    int i,j,k,matsize;

    printf("Enter the size of the matrix(i.e. value of 'n' as size is
nXn):");
    scanf("%d",&matsize);

    A=(float **)malloc(matsize*sizeof(float *));            //allocate memory
dynamically for matrix A(matsize X matsize)
    for(i=0;i<matsize;i++)
        A[i]=(float *)malloc(matsize*sizeof(float));

    I=(float **)malloc(matsize*sizeof(float *));            //memory allocation for
indentity matrix I(matsize X matsize)
    for(i=0;i<matsize;i++)
        I[i]=(float *)malloc(matsize*sizeof(float));

    printf("Enter the matrix: ");                           // ask the user for matrix A
    for(i=0;i<matsize;i++)
        for(j=0;j<matsize;j++)
            scanf("%f",&A[i][j]);

    for(i=0;i<matsize;i++)                                  //automatically initialize the unit matrix, e.g.
        for(j=0;j<matsize;j++)                              //  -       -
            if(i==j)                                        // | 1  0  0 |
                I[i][j]=1;                                  // | 0  1  0 |
            else                                            // | 0  0  1 |
                I[i][j]=0;                                  //  -       -
/*---------------LoGiC starts here------------------*/      //procedure // to make the matrix A to unit matrix

    for(k=0;k<matsize;k++)                                  //by some row operations,and the same row operations of
    {                                                       //Unit mat. I gives the inverse of matrix A
        temp=A[k][k];                   //'temp'  
        // stores the A[k][k] value so that A[k][k]  will not change
        for(j=0;j<matsize;j++)      //during the operation //A[i] //[j]/=A[k][k]  when i=j=k
        {
            A[k][j]/=temp;                                  //it performs // the following row operations to make A to unit matrix
            I[k][j]/=temp;                                  //R0=R0/A[0][0],similarly for I also
R0=R0/A[0][0]
        }                                                   //R1=R1-R0*A[1][0] similarly for I
        for(i=0;i<matsize;i++)                              //R2=R2-R0*A[2][0]      ,,
        {
            temp=A[i][k];                       //R1=R1/A[1][1]
            for(j=0;j<matsize;j++)             //R0=R0-R1*A[0][1]
            {                                   //R2=R2-R1*A[2][1]
                if(i==k)
                    break;                      //R2=R2/A[2][2]
                A[i][j]-=A[k][j]*temp;          //R0=R0-R2*A[0][2]
                I[i][j]-=I[k][j]*temp;          //R1=R1-R2*A[1][2]
            }
        }
    }
/*---------------LoGiC ends here--------------------*/
    printf("The inverse of the matrix is: ");               //Print the //matrix I that now contains the inverse of mat. A
    for(i=0;i<matsize;i++)
    {
        for(j=0;j<matsize;j++)
            printf("%f  ",I[i][j]);
        printf(" ");
    }
    return 0;
}

Answer 2

对角线元素必须首先缩放到1，然后才能将下三角形元素归零(第二个嵌套循环)。结果是对角线包含0，=>，不存在逆，或者我们得到一个行梯形。通过对角线上的反向迭代，我们可以减少它并得到逆。代码远非最优。当你有FPU的时候，对于这样的数值计算来说，double可能是比float更好的选择。

请注意，零化子矩阵、行交换等可能会被更优化的解决方案所取代。Matrix是一个自定义类型，IsFloat0是一个自定义函数，但都应该清楚命名和上下文。享受代码：

const uint sz = 4;
Matrix< double > mx;
mx.Resize( 2 * sz, sz );
mx.Zero();
for( uint rdx = 0; rdx < mx.NumRow(); ++rdx )
{
        mx( rdx, rdx + mx.NumRow() ) = 1.0; // eye
}
mx( 0, 0 ) = 1.0; mx( 0, 1 ) = 1.0; mx( 0, 2 ) = 1.0; mx( 0, 3 ) = 0.0;
mx( 1, 0 ) = 1.0; mx( 1, 1 ) = 1.0; mx( 1, 2 ) = 2.0; mx( 1, 3 ) = 0.0;
mx( 2, 0 ) = 1.0; mx( 2, 1 ) = 2.0; mx( 2, 2 ) = 0.0; mx( 2, 3 ) = 1.0;
mx( 3, 0 ) = 1.0; mx( 3, 1 ) = 2.0; mx( 3, 2 ) = 0.0; mx( 3, 3 ) = 2.0;

// pivot iteration
uint idx;
for( idx = 0; idx < sz; ++idx )
{
    // search for non-0 pivot
    uint sdx = sz;
    for( uint rdx = idx; rdx < sz; ++rdx )
    {
        if( !Util::IsFloat0( mx( rdx, idx ) ) ) { sdx = rdx; rdx = sz - 1; }
    }
    if( sdx < sz )
    {
        // swap rows
        if( idx != sdx )
        {
            for( uint cdx = 0; cdx < ( sz << 1 ); ++cdx )
            {
                double swp;
                swp = mx( idx, cdx );
                mx( idx, cdx ) = mx( sdx, cdx );
                mx( sdx, cdx ) = swp;
            }
        }
        // 1 pivot and 0 col
        {
            double sc = 1.0 / mx( idx, idx );
            for( uint cdx = 0; cdx < ( sz << 1 ); ++cdx )
            {
                mx( idx, cdx ) *= sc; // 1
            }

            for( uint rdx = 1 + idx; rdx < sz; ++rdx )
            {
                double sd = mx( rdx, idx );
                for( uint cdx = 0; cdx < ( sz << 1 ); ++cdx )
                {
                    mx( rdx, cdx ) -= sd * mx( idx, cdx ); // 0
                }
            }
        }
    }
    else { idx = sz; }
}
if( sz < idx ) { mx.Zero(); }
else
{
    for( idx = 0; idx < sz; ++idx )
    {
        uint ydx = sz - 1 - idx;
        for( uint rdx = 0; rdx < ydx; ++rdx )
        {
            double sc = mx( rdx, ydx );
            for( uint cdx = 0; cdx < ( sz << 1 ); ++cdx )
            {
                mx( rdx, cdx ) -= sc * mx( ydx, cdx ); // 0
            }
        }
    }
}