CODEWITHSUNDEEP
c++intel-mkleigen3
Tom de Geus
Tom de Geus

Reputation: 5965

Use Eigen + Intel MKL + Pardiso

I am trying to use Eigen's support of MKL and Pardiso (see example below). I have used the Intel link line advisor to come up with the compiler options but everything I'm trying is unsuccessful. In particular, compiling with:

$ icpc -I${HOME}/src/eigen -DEIGEN_USE_MKL_ALL -DMKL_ILP64 -I${MKLROOT}/include main.cpp -L${MKLROOT}/lib/intel64 -lmkl_intel_ilp64 -lmkl_sequential -lmkl_core -lpthread -lm -ldl

Results in the following type error:

~/src/eigen/Eigen/src/PardisoSupport/PardisoSupport.h(50): error: argument of type "int *" is incompatible with parameter of type "const long long *" ::pardiso(pt, &maxfct, &mnum, &type, &phase, &n, a, ia, ja, perm, &nrhs, iparm, &msglvl, b, x, &error);

While compiling without the 64-bit support, using:

$ icpc -I${HOME}/src/eigen -DEIGEN_USE_MKL_ALL -I${MKLROOT}/include main.cpp -L${MKLROOT}/lib/intel64 -lmkl_sequential -lmkl_core -lpthread -lm -ldl

Results in:

main.cpp:(.text+0x759): undefined reference to `pardiso'

I have also tried, equally unsuccessfully, to compile with OpenMP.

How should I compile?

(I am on Intel 2016.3.210 and the latest Eigen, 3.3.3)

Example (main.cpp)

This example solves Poisson's equation, with two Dirichlet boundary conditions (u_{0} = u_{N+1} = 0)

#include <iostream>
#include <iomanip>
#include <math.h>
#include <Eigen/Dense>
#include <Eigen/PardisoSupport>

int main()
{
  typedef Eigen::SparseMatrix<double> SpMat;
  typedef Eigen::Triplet     <double> Trip;
  typedef Eigen::PardisoLU   <SpMat > Solver;

  size_t N = 11;

  SpMat A(N,N);
  Eigen::VectorXd f(N);

  f         *= 0.0;
  f((N-1)/2) = 1.0;

  std::vector<Trip> Atr;

  Atr.push_back(Trip(0,0,+2.0));
  Atr.push_back(Trip(0,1,-1.0));

  for ( size_t i=1; i<N-1; ++i ) {
    Atr.push_back(Trip(i,i-1,-1.0));
    Atr.push_back(Trip(i,i  ,+2.0));
    Atr.push_back(Trip(i,i+1,-1.0));
  }

  Atr.push_back(Trip(N-1,N-2,-1.0));
  Atr.push_back(Trip(N-1,N-1,+2.0));

  A.setFromTriplets(Atr.begin(),Atr.end());

  Solver solver;
  solver.analyzePattern(A);
  solver.factorize(A);

  Eigen::VectorXd u = solver.solve(f);

  return 0;
}

Upvotes: 0

Views: 2286

Answers (1)

Tom de Geus
Tom de Geus

Reputation: 5965

It turns out the the eigen3 documentation is very clear about this:

Using Intel MKL through Eigen is easy:

...

  1. on a 64bits system, you must use the LP64 interface (not the ILP64 one)

I have now successfully compiled with 

icpc -I${HOME}/src/eigen -DEIGEN_USE_MKL_ALL -DMKL_LP64 -I${MKLROOT}/include main.cpp -L${MKLROOT}/lib/intel64 -lmkl_intel_lp64 -lmkl_sequential -lmkl_core -lpthread -lm -ldl

Upvotes: 1

Related Questions