使用具有64位整数的稀疏矩阵作为StorageIndex的Eigen Pardiso

Question

我最近在使用Eigen Pardiso开发代码时遇到了一些问题。

以前，我使用Eigen Pardiso来求解二维泊松方程。矩阵大小可以用C/C++32位整数来描述。我使用的代码是这样的：

#define EIGEN_USE_MKL_ALL
typedef Eigen::SparseMatrix<double, Eigen::ColMajor> SpMat;
typedef Eigen::Triplet<double> Trp;
......
int nx;
int ny;
SpMat A(nx*ny, nx*ny);
Eigen::VectorXd b=Eigen::VectorXd::Zeros(nx*ny);
Eigen::VectorXd sol=Eigen::VectorXd::Zeros(nx*ny);
Eigen::PardisoLU<SpMat> *ptrLUsolver;
ptrLUsolver = new Eigen::PardisoLU<SpMat>();
// Set-up coefficients for A and set RHS vector b
......
//
(*ptrLUsolver).compute(A);
sol = (*ptrLUsolver).solve(b);
if (ptrLUsolver != NULL)
{
delete ptrLUsolver;
}
......

这对我来说很好。

然而，我最近需要修改我的代码来处理更大的3D矩阵。上述代码中的nx*ny需要更改为nx*ny*nz，并且该值大于C/C++中的32位int。因此，这些变量需要在C/C++中转换为64位的CCD_ 4。上面代码的修改代码如下所示：

#define EIGEN_USE_MKL_ALL
typedef Eigen::SparseMatrix<double, Eigen::ColMajor, int64_t> SpMat;
typedef Eigen::Triplet<double, int64_t> Trp;
......
int64_t nx;
int64_t ny;
int64_t nz;
SpMat A(nx*ny*nz, nx*ny*nz);
Eigen::VectorXd b=Eigen::VectorXd::Zeros(nx*ny*nz);
Eigen::VectorXd sol=Eigen::VectorXd::Zeros(nx*ny*nz);
Eigen::PardisoLU<SpMat> *ptrLUsolver;
ptrLUsolver = new Eigen::PardisoLU<SpMat>();
// Set-up coefficients for A and set RHS vector b
......
//
(*ptrLUsolver).compute(A);
sol = (*ptrLUsolver).solve(b);
if (ptrLUsolver != NULL)
{
delete ptrLUsolver;
}

现在问题来了。似乎通过将SparseMatrix类中的默认StorageIndex更改为int64_t，代码根本无法编译。错误信息如下：

error: cannot convert ‘long int*’ to ‘const int*’
::pardiso(pt, &maxfct, &mnum, &type, &phase, &n, a, ia, ja, perm, &nrhs, iparm, &msglvl, b, x, &error);
~~~~~~~~~^~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

然后我来谈谈Eigen的源代码。似乎有这样的东西：

namespace internal
{
template<typename IndexType>
struct pardiso_run_selector
{
static IndexType run( _MKL_DSS_HANDLE_t pt, IndexType maxfct, IndexType mnum, IndexType type, IndexType phase, IndexType n, void *a,
IndexType *ia, IndexType *ja, IndexType *perm, IndexType nrhs, IndexType *iparm, IndexType msglvl, void *b, void *x)
{
IndexType error = 0;
::pardiso(pt, &maxfct, &mnum, &type, &phase, &n, a, ia, ja, perm, &nrhs, iparm, &msglvl, b, x, &error);
return error;
}
};
template<>
struct pardiso_run_selector<long long int>
{
typedef long long int IndexType;
static IndexType run( _MKL_DSS_HANDLE_t pt, IndexType maxfct, IndexType mnum, IndexType type, IndexType phase, IndexType n, void *a,
IndexType *ia, IndexType *ja, IndexType *perm, IndexType nrhs, IndexType *iparm, IndexType msglvl, void *b, void *x)
{
IndexType error = 0;
::pardiso_64(pt, &maxfct, &mnum, &type, &phase, &n, a, ia, ja, perm, &nrhs, iparm, &msglvl, b, x, &error);
return error;
}
};

根据错误信息，我认为代码调用了pardiso函数。所以我认为代码可能需要找到一些方法来使用pardiso_64函数。我正试图在上面所示代码的上述修改版本上将int64_t更改为long long int。然后编译器显示错误：

eigen-3.3.9/Eigen/src/Core/util/XprHelper.h:833:24: error: static assertion failed: YOU_MIXED_DIFFERENT_NUMERIC_TYPES__YOU_NEED_TO_USE_THE_CAST_METHOD_OF_MATRIXBASE_TO_CAST_NUMERIC_TYPES_EXPLICITLY
EIGEN_STATIC_ASSERT((Eigen::internal::has_ReturnType<ScalarBinaryOpTraits<LHS, RHS,BINOP> >::value), 
eigen-3.3.9/Eigen/src/Core/util/StaticAssert.h:33:54: note: in definition of macro ‘EIGEN_STATIC_ASSERT’
#define EIGEN_STATIC_ASSERT(X,MSG) static_assert(X,#MSG);
eigen-3.3.9/Eigen/src/Core/AssignEvaluator.h:834:3: note: in expansion of macro ‘EIGEN_CHECK_BINARY_COMPATIBILIY’
EIGEN_CHECK_BINARY_COMPATIBILIY(Func,typename ActualDstTypeCleaned::Scalar,typename Src::Scalar);
^~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
eigen-3.3.9/Eigen/src/SparseCore/SparseSelfAdjointView.h:641:59: error: cannot bind non-const lvalue reference of type ‘Eigen::SparseMatrix<double, 0, long long int>&’ to an rvalue of type ‘Eigen::SparseMatrix<double, 0, long long int>’
internal::permute_symm_to_fullsymm<Mode>(src.matrix(),tmp,src.perm().indices().data());

我可以问一下如何将特征中的Pardiso与int64_t一起使用吗？提前谢谢。

我使用的特征是Eigen 3.3.9。我使用的编译器是GCC 8.4.0。

Answer 1

我认为在64位系统上，Eigen使用long int作为索引和大小的默认类型。但是pardiso_64需要每个参数都和int一样长

https://software.intel.com/content/www/us/en/develop/documentation/onemkl-developer-reference-c/top/sparse-solver-routines/onemkl-pardiso-parallel-direct-sparse-solver-interface/pardiso-64.html

在Eigen文档中，默认情况下，Eigen_DEFAULT_DENSE_INDEX_TYPE设置为std:：ptrdiff_t。这个std:：ptrdiff_t用于指针运算和数组索引。使用其他类型(如int(的程序可能会在64位系统上失败，例如，当索引超过int_MAX或依赖于32位模块算术时。

因此，请尝试将EIGEN_DEFAULT_DENSE_INDEX_TYPE更改为long-long-int.

https://eigen.tuxfamily.org/dox/TopicPreprocessorDirectives.html