1

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

以前，我使用 Eigen Pardiso 来求解 2D Poisson 方程。矩阵大小可以用 C/C++ 32 位 int 数来描述。我使用的代码是这样的：

#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，这个值比intC/C++ 中的 32-bit 大。int64_t所以这些变量需要在 C/C++ 中改为 64 位。列出了上面代码的修改代码：

#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, &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());

我可以问如何使用 Eigen 中的 Pardisoint64_t吗？提前致谢。

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

4

1 回答 1

0

我认为在 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-接口/pardiso-64.html

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

所以尝试将 EIGEN_DEFAULT_DENSE_INDEX_TYPE 更改为 long long int。

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

于 2021-07-01T05:19:20.773 回答