如何解决使用带有 64 位整数作为 StorageIndex 的 SparseMatrix 的 Eigen Pardiso
我最近在使用 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 位 int64_t。上面代码修改后的代码如下：
#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;
}
现在问题来了。好像把StorageIndex类的默认SparseMatrix改成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,&error);
return error;
}
};
template<>
struct pardiso_run_selector<long long int>
{
typedef long long int IndexType;
static IndexType run( _MKL_DSS_HANDLE_t pt,void *x)
{
IndexType error = 0;
::pardiso_64(pt,&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,long long int>&’ to an rvalue of type ‘Eigen::SparseMatrix<double,long long int>’
internal::permute_symm_to_fullsymm<Mode>(src.matrix(),tmp,src.perm().indices().data());
请问如何在 Eigen 中使用 Pardiso 和 int64_t？提前致谢。
我使用的特征是 Eigen 3.3.9。
我使用的编译器是 GCC 8.4.0。
解决方法
我认为在 64 位系统上，Eigen 使用 long int 作为索引和大小的默认类型。但是 pardiso_64 需要每个参数都和 long long int 一样长。
在 Eigen 文档中，EIGEN_DEFAULT_DENSE_INDEX_TYPE 默认设置为 std::ptrdiff_t。这个 std::ptrdiff_t 用于指针运算和数组索引。使用其他类型（例如 int）的程序可能会失败，例如当索引超过 INT_MAX 或者它依赖于 32 位模运算时，64 位系统。
因此尝试将 EIGEN_DEFAULT_DENSE_INDEX_TYPE 更改为 long long int。
https://eigen.tuxfamily.org/dox/TopicPreprocessorDirectives.html
版权声明：本文内容由互联网用户自发贡献，该文观点与技术仅代表作者本人。本站仅提供信息存储空间服务，不拥有所有权，不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容， 请发送邮件至 dio@foxmail.com 举报，一经查实，本站将立刻删除。
