这将是一个非常普遍的问题。

我在 Python 中有一个 3D 有限元代码，我想扩展它以处理“大”问题（全球系统中约 10^8 个未知数）。现在我正在使用该scipy.sparse库，它为迭代求解器提供了不错的性能，但我发现了以下问题：

• 对于大于 10^6 个未知数的问题，我很快就会遇到内存限制
• 我正在解决的线性系统是对称正定的，但scipy.sparse似乎没有一种了解对称性的存储格式，所以我认为我存储的条目比必要的多得多。
• 甚至在求解全局系统之前，计算元素局部贡献（保存在多维numpy数组中）就已经接近内存需求。

因此，我似乎很清楚我需要编写自己的无矩阵迭代求解器，或者使用具有该功能的外部库。我的问题：

• 我知道它scipy.sparse包含较低级别的例程，例如 lapack。我需要写的基本上是重载稀疏矩阵向量乘法。如果我用 Python（或者可能用 C 语言并用 Python 包装）编写它，我是否有希望获得不错的性能，或者只是需要能够调用这些较低级别的库？我对此没有直觉，也不想花 3 周时间编写我自己的太慢的无矩阵 A*x 例程。
• 用 Python 编写这个例程是个好主意吗？还是我需要用 C 编写并包装它？
• 是否存在支持此功能的库？这似乎不太可能，因为有必要了解如何在不组装 A 的情况下计算 A*x 的代码特定知识。
• 我可以编写一个核外求解器，而不是编写一个无矩阵的例程吗？这种方法是否适用于像我这样的问题？
• 我可以访问具有许多多核计算节点的集群。最终我想并行化这个实现。实施指南是否有好的工具或参考资料？

或者有没有一种很好的方法来处理上述问题scipy？例如，是否可以提供scipy.sparse.linalg.cg一个指向计算 A * x 的函数的指针（显然 cg 例程受 A * x 速度的限制，但最好避免自己编写 CG、GMRES 等）？