使用具有64位整数作为StorageIndex的SparseMatrix的Eigen Pardiso

Question

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

之前，我使用Eigen Pardiso来求解2D Poisson方程。矩阵大小可以用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;
}

现在问题来了。似乎通过将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());

我可以问一下如何在int64_t中使用Eigen中的Pardiso吗？提前谢谢。

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

Answer 1

我认为在64位系统上，Eigen使用long int作为索引和大小的默认类型。但是pardiso_64需要每个参数都是long long 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_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