2x2矩阵的数值稳定逆

Question

在我用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来确定我是否可以avoid the inversion。