相关疑难解决方法(0)

在我在C中工作的数值解算器中,我需要反转2x2矩阵,然后在右侧乘以另一个矩阵:

C = B . inv(A)

我一直在使用倒置2x2矩阵的以下定义:

a = A[0][0];
b = A[0][1];
c = A[1][0];
d = A[1][1];
invA[0][0] = d/(a*d-b*c);
invA[0][1] = -b/(a*d-b*c);
invA[1][0] = -c/(a*d-b*c);
invA[1][1] = a/(a*d-b*c);

在我解算器的前几次迭代中,这似乎给出了正确的答案,然而,经过几个步骤后,事情开始增长并最终爆炸.

现在,与使用SciPy的实现相比,我发现相同的数学不会爆炸.我能找到的唯一区别是SciPy代码使用scipy.linalg.inv(),它在内部使用LAPACK来执行反转.

当我用inv()上面的计算替换调用时,Python版本确实会爆炸,所以我很确定这是问题所在.计算中的微小差异正在蔓延,这使我相信它是一个数值问题 - 对于反演操作并不完全令人惊讶.

我正在使用双精度浮点数(64位),希望数值问题不会成为问题,但显然情况并非如此.

但是:我想在我的C代码中解决这个问题,而不需要调用像LAPACK这样的库,因为将它移植到纯C的全部原因是让它在目标系统上运行.此外,我想了解这个问题,而不仅仅是呼唤黑匣子.最终,如果可能的话,我也希望它以单精度运行.

所以,我的问题是,对于这样一个小矩阵,有一个数值上更稳定的方法来计算A的倒数吗？

谢谢.

编辑:目前试图弄清楚我是否可以通过解决来避免反演C.