c++ - 并行场景中的 Eigen::PardisoLU (MKL)

Question

我正在尝试将 Eigen（3.3.1 或 3.3.2）与 INTEL 的 MKL PardisoLU 求解器一起使用。我正在使用 EIGEN_USE_MKL_ALL 标志并安装了 IntelSWTools 2017。目前仅在带有 VS C++ 的 Win10 上尝试过，但这不是目标平台。

求解器本身看起来比 SparseLU 工作得快得多，这很好（我对使用非常大的矩阵很感兴趣）。但是在并行场景中，它不能正确解决我的问题。SparseLU 没有这样的问题，或者当我删除上层并行性时（评论pragma parallel for下面）。

PardisoLU 本身应该是线程安全的，那我做错了什么？

该场景非常简短，如下所示

#include <Eigen/PardisoSupport>

#define N = 100000                   

typedef Eigen::Matrix<Complex, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor> FullMatrix;

typedef Eigen::SparseMatrix<Complex, Eigen::ColMajor, long long> SparseMatrix; 

typedef Eigen::PardisoLU<SparseMatrix> Solver;
//typedef Eigen::SparseLU<SparseMatrix, Eigen::COLAMDOrdering< long long>> Solver;

Eigen::initParallel();

SparseMatrix A = CreateA();

Solver S;
S.analyzePattern(A);
S.factorize(A);

FullMatrix Q(A.rows(),N);

#pragma omp parallel for  
for (int n = 0; n < N; ++n)
{
     SparseMatrix Rhs = getRhs(n);
     Q.col(n) = S.Solve(Rhs);
}