用C方阵反转 [英] Square Matrix Inversion in C

问题描述

我写了一个反转功能为n * n的方阵。

I wrote a inversion function for an n*n square matrix.

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

这正常工作在几乎所有情况下，但在某些情况下，像此处显示的失败：

This works fine in almost all cases but is failing in some cases like the ones shown here:

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

还有什么我俯瞰的情况？

What could be the case I am overlooking?

谢谢！

推荐答案

对角线elems必须先归零下三角elems（你的第二个嵌套循环）之前缩放到1秒。事实证明，对角线包含0 =>不存在逆或者我们得到一个行梯形式。随着对角线逆转迭代，我们可以减少它，并得到逆。
在code是很不理想。当你有FPU，双能比浮动此类数值计算一个更好的选择。

Diagonal elems have to be scaled to 1s first before zeroing lower triangle elems (your second nested loop). It turns out that diagonal contains 0 => no inverse exists OR we get a row echelon form. With a reversed iteration on diagonal we can reduce it and get the inverse. The code is far from optimal. When you have FPU, double could be a better choice than float for such numerical computations.

需要注意的是零子矩阵，排掉期等可以通过更优化的解决方案所取代。 Matrix是一个自定义类型和IsFloat0是一个自定义的功能，但应该明确的命名和背景。享受code：

Note that zeroing submatrices, row swaps etc. could be replaced by far more optimal solutions. Matrix is a custom type and IsFloat0 is a custom function but all should be clear of naming and context. Enjoy code:

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

这篇关于用C方阵反转的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！