C中的方阵反演

我为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==(jn)) 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 。

使用Gauss Jordan算法

 #include #include int main() { float **A,**I,temp; int i,j,k,matsize; printf("Enter the size of the matrix(ie 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 

Answer 2

在将下三角形元素（第二个嵌套循环）归零之前，必须首先将对角元素缩放为1。 事实certificate，对角线包含0 =>没有逆存在或者我们得到行梯形forms。 通过对角线上的反向迭代，我们可以减少它并得到逆。 代码远非最佳。 当你有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 } } } }