Using Eigen Pardiso with SparseMatrix which has 64 bit integer as StorageIndex

Question

I recently got some problems when developing my code using Eigen Pardiso.

Previously, I was using Eigen Pardiso to solve the 2D Poisson equation. The Matrix size is OK to be described by the C/C++ 32 bit int number. The code I was using is like this:

#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;
}
......

This works well for me.

However, I recently need to modify my code to work on a larger 3D matrix. The nx*ny in the above code needs to be changed into nx*ny*nz and this value is larger than the 32-bit int in the C/C++. So those variables need to be changed into 64 bit int64_t in the C/C++. The modified code from the code above is listed:

#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;
}

Now the problem comes. It seems like by changing the default StorageIndex in the SparseMatrix class into int64_t, the code can not be compiled at all. The error information is like this:

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);
       ~~~~~~~~~^~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

Then I come to the source code of Eigen. It seems like there is something like:

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;
    }
  };

From the error information, I think the code calls the pardiso function. So I think probably the code needs to find some way to use the pardiso_64 function. I am trying to change the int64_t into long long int on the above-modified version of the code shown above. Then the compiler shows errors:

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());

Can I ask how to use the Pardiso in Eigen with int64_t? Thanks in advance.

The Eigen I am using is Eigen 3.3.9. The compiler I am using is GCC 8.4.0.

@김선달 Hi, first thanks for the comment. I think the results for the `std::is_same::value` is false. — Chen Cui, Jun 23 '21 at 05:26
Then that must be a problem. Change your `SpMat` and `TpMat` typedefs to use `long long int` instead of `int64_t` and give it a try. — 김선달, Jun 23 '21 at 05:31
@김선달 Thanks for the comment. I tried to substitute the `int64_t` to `long long int` but the compiler still shows the same as shown in the last part of the initial post. — Chen Cui, Jun 23 '21 at 05:56

score 0 · Answer 1 · answered Jul 01 '21 at 05:19

I think on 64 bits system Eigen uses long int as the default type for indexes and sizes. But the pardiso_64 needs every argument as 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

In the Eigen documentation, EIGEN_DEFAULT_DENSE_INDEX_TYPE is set to std::ptrdiff_t by default. And this std::ptrdiff_t is used for pointer arithmetic and array indexing. Programs that use other types, such as int, may fail on, e.g. 64-bit systems when the index exceeds INT_MAX or if it relies on 32-bit modular arithmetic.

So try changing the EIGEN_DEFAULT_DENSE_INDEX_TYPE to long long int.

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

Thanks a lot! I will give that a try! – Chen Cui Jul 01 '21 at 23:12 — Chen Cui, Jul 01 '21 at 23:12

Using Eigen Pardiso with SparseMatrix which has 64 bit integer as StorageIndex

1 Answers1