boost::numeric::ublas::matrix LU求逆方法可靠性咨询

需要对3030左右规模的双精度浮点数方阵求逆。
最初计划自行用C++实现LU分解方法，后续了解到boost::numeric::ublas::matrix库，为避免重复开发，调用该库的相关函数，按顺序执行lu_factorize()和lu_substitute()完成矩阵求逆操作。
手动验证该求逆逻辑（仅调用上述两个boost函数）的可靠性时，3阶、4阶简单方阵的测试结果均符合预期。
但对3030矩阵A求逆得到A⁻¹后，计算乘积A*A⁻¹得到的矩阵仅对角线元素为1，其余位置均为极小值，结果片段如下：

1         -1.5e-16  -5.1e-20  2.4e-19  
 0         1         0         -5.4e-20  
 1.1e-16   1.1e-16   1         -1.4e-19  
 6.9e-17   0         -3.3e-17  1         

无法判断这些非对角线位置的极小值是矩阵库近似计算0产生的误差，还是LU分解过程引入的误差。

• 是否有开发者使用boost ublas时遇到过同类问题？
• 该库目前是否可靠、是否已过时？
• 是否可以获取这些算法的源代码？

开发环境：C++11、gcc/8.2.0、boost/1.76.0

你看到的非对角线极小值是双精度浮点数计算的正常数值误差，不是boost ublas的实现bug。
双精度浮点数的有效精度约为15_{17位十进制数，机器epsilon（即1和下一个可表示双精度数的差值）约为2.2e-16。30阶矩阵的LU分解、求逆、矩阵乘法全流程会累积上百次浮点运算，最终非对角元素出现1e-16}1e-19量级的偏差完全在预期范围内，这个量级的数值和理论值0的差异已经小于双精度能分辨的最小误差，完全可以判定为数值计算下的等价0，不需要额外处理。如果业务场景对精度要求更高，可以考虑切换到长双精度类型，不过常规科学计算场景下这个精度完全够用。

针对三个具体问题的说明：

• 这类浮点数误差是所有线性代数计算库的共性问题，不是ublas独有，用任何库做30阶双精度矩阵求逆都会出现同量级的偏差，很多开发者第一次接触浮点数矩阵计算都会遇到同样的疑惑。
• ublas是Boost官方维护的线性代数组件，基础LU分解、回代求逆的逻辑经过了二十多年的生产环境验证，不存在可靠性问题，你在3、4阶矩阵上的测试结果也能佐证核心逻辑正确。它确实接口设计偏老旧，没有针对现代CPU的SIMD、多线程做深度优化，性能比Eigen、OpenBLAS这类专门优化过的线性代数库差不少，但不存在“过时到无法使用”的问题，你用到的30阶小规模矩阵场景，性能差异根本感知不到。
• 这些算法的源代码完全公开，不需要额外下载，你本地安装的boost 1.76.0目录下就能找到对应实现：lu_factorize和lu_substitute的全部实现在boost/numeric/ublas/lu.hpp头文件中，直接打开就能阅读完整源码，ublas的大部分实现都是头文件内联的，没有编译成闭源二进制。

补充一个实用的验证方法：不用逐元素核对非对角项，直接计算A*A⁻¹和单位矩阵的Frobenius范数差即可，30阶双精度矩阵下这个差值在1e-14量级以内就属于完全正常的结果。

内容的提问来源于stack exchange，提问作者SoleGoodman