鉴于 UMFPACK 内存限制，如何解决大型系统？

计算科学 线性代数 C++ 布拉斯 内存管理 umfpack
2021-12-22 22:32:33

我正在尝试使用 UMFPACK 和与 C++ 的增强数值绑定来求解 3D 热扩散的方程组（ A x = b）（即每个方程最多有 7 个项，不包括常数“b”项）。

例如：

namespace umf = boost::numeric::bindings::umfpack;

umf::umf_solve(A,x,b);

我可以使用大约 200MB 的内存为 200K+ 组方程填充A矩阵和“b”向量。然后在 UMFPACK 例程中，我的程序使用的内存跃升至 2.3GB+。

此外，在约 211K 方程以上的某个点，例程将“x”的所有零传回，而不会引发错误。我认为这与某种内存限制有关。

这在我测试过的每台机器上都是一样的（我的 MacBook、我的 Ubuntu VM 和 Red Hat 超级计算机上的一个节点）。

我的问题是：

  1. 为什么 UMFPACK 返回零而没有错误？
  2. 为什么 UMFPACK 需要这么多内存？
  3. 我有哪些选项可以为更大的方程组求解系统？
2个回答

我不能给你一个好的答案 1，但我可以给你体面的答案 2 和 3。

  1. 我对 UMFPACK 的 Boost 接口不是很熟悉。在 C 接口中，通常调用将分配内存的 UMFPACK 例程。如果由于没有足够的可用内存而无法分配内存，UMFPACK 将返回一个空指针。然后，您可以测试返回空指针以处理错误的例程；这种方法在 C 编程中是标准的。查看文档，您似乎必须检查每个函数的返回值；看起来图书馆没有抛出任何异常。如果您不检查返回值，而仅检查 UMFPACK 返回的解决方案，则可以解释您观察到的行为。

  2. UMFPACK 是一个稀疏的直接求解器；矩阵执行稀疏 LU 分解。随着变大，更多 LU 因式分解的条目在结构上变得非零——它们被“填充”，这种现象称为填充——结构上的非零条目在内存中显式表示，因为它们可能在数值上非零. 填充是你内存不足的原因；内存限制往往使稀疏的直接求解器对超过大约 100K 未知数的问题没有吸引力。AA

  3. 使用迭代法求解稀疏线性方程，而不是稀疏直接法。这些方法不一定要求您显式存储；它们通常只需要一种计算任何的方法。为了获得最佳结果，您应该为迭代方法提供一个预处理器。对于扩散问题，多重网格是一个很好的预处理器；多重网格方法也值得在没有扩散方程预处理器的情况下使用。其他人可能会就用于此类问题的迭代方法提供更好的建议。AAxx

Geoff 提到了稀疏直接求解器的局限性。要解决大小为的 3D 问题，最佳直接求解器将需要内存和时间。有关更多详细信息，请参阅此问题。N=n3N4/3N2

但听起来你有一个结构化的网格扩散问题。如果扩散系数是平滑的（或恒定的），您可以使用著名的 Full Multigrid 算法（由 Achi Brandt 在 1970 年代开发）使用大约四个浮点值的内存和每个网格点 50 次浮点数来解决问题。如果扩散系数非常粗糙，最容易使用代数多重网格方法。根据您的需要，PETSc和其他软件包中提供了几何和代数多重网格方法。[免责声明：我是 PETSc 的开发人员。]

其它你可能感兴趣的问题
上一篇大矩阵上的Matlab预处理共轭梯度 下一篇如何判断哪个（顺序或并行）版本的英特尔 MKL 已链接？