使用 eigen/Intel MKL 求解稀疏线性系统答案

【问题标题】：Solving sparse linear system with eigen/Intel MKL使用 eigen/Intel MKL 求解稀疏线性系统
【发布时间】：2023-03-13 07:54:01
【问题描述】：

我想通过使用 C++ 中的 Eigen 库来求解 Ax=b 方程，其中 A 是稀疏矩阵 (1,964,568 x 1,964,568 nnz=75256446)，b 也是稀疏矩阵 (1,964,568 x 1,964,568 nnz= 25354926)。

起初我试图使用 Eigen Sparse LU 来解决我的问题，但几个小时后我的内存不足（我有 128GB RAM）。在此之后，我将 INTEL MKL 库与 Pardiso 求解器一起包含在内。即使这样，我也无法解决我的问题。也许有人有一些技巧可以解决我的问题？

#define EIGEN_USE_MKL_ALL

#include <Eigen/Sparse>
#include <unsupported/Eigen/SparseExtra>
#include <iostream>
#include <Eigen/OrderingMethods>
#include <Eigen/PardisoSupport>


typedef Eigen::SparseMatrix<double>SpMat; 
typedef Eigen::COLAMDOrdering<int>Order;

int main()
{

    SpMat A;
    SpMat B;
    SpMat X;

    Eigen::loadMarket(A, "MatK.mtx");
    Eigen::loadMarket(B, "MatM.txt");

    A.makeCompressed();
    B.makeCompressed();

    Eigen::PardisoLU<SpMat>solver;

    solver.analyzePattern(A);

    solver.factorize(A);

    X = solver.solve(B);


}

我可以编译我的代码并运行它。我只需要更好的性能和更少的内存。

【问题讨论】：

PARDISO 是一个直接求解器，可能需要大量额外内存。尝试使用一些 iterative 求解器。这些可能还需要一些内存，例如，用于预处理和存储 Krylov 子空间向量，但通常不如直接求解器那么多。 Eigen 和英特尔 MKL 都包含一些迭代求解器。顺便说一句，您使用 64 位二进制文件和库吗？
所以您正在尝试解决A * X = B 的X 问题，其中A 和B 都是巨大的稀疏矩阵？你假设解决方案X 也是一个稀疏矩阵？通常的情况是x 和b 是向量。不确定是否可以使用通常的求解器获得稀疏解矩阵X。我可能宁愿解决B 的每一列b_i 并将结果x_i 组合成一个矩阵X。
是的，我注意到我必须对 B 矩阵进行切片。我想要解决的问题是：X=inv(A)*B。因为我不想计算逆我把它写成 AX=B 并解决这个问题。
填充量（以及所需的内存量）很大程度上取决于矩阵的结构（以及似乎无法为Eigen::ParadisoLU 配置的排序方法） .您的矩阵中是否有一些可以利用的特殊结构？要获得明确的答案，您不仅应该提供来源，还应该提供重现问题所需的所有数据（请参阅minimal reproducible example）。
我看不出有任何理由假设 X 是稀疏的，即使是这样，我也怀疑求解器会返回 SparseMatrix。

标签： c++ sparse-matrix eigen pardiso

【解决方案1】：

因为我还没有找到任何方法来解决这个问题。我试图改变我的 RHS 矩阵，因为大量的列是一个大问题。我在我的算法中看到我的 RHS 乘以另一个矩阵，这将列减少到 10。有了这个，就可以使用英特尔 Pardiso LDLT 求解器而不是 LU 来求解。迭代求解器需要太多迭代才能收敛。

【讨论】：