Hello, I'm building my application using VS2015, Eigen 3.3.4, with next flags

EIGEN_USE_LAPACKE
EIGEN_USE_BLAS

First,
I have changed in lapacke.h
#define lapack_complex_float float _Complex
#define lapack_complex_double double _Complex

to

#define lapack_complex_float _Fcomplex
#define lapack_complex_double _Dcomplex

otherwise I'll always obtain redifinition errors for _Complex type

Second,
I can link only using Intel MKL libraries.

Linking with original LAPACK.lib, BLAS.lib & LAPACKE.lib returns an error

error LNK2019: unresolved external symbol dgemm_ referenced in function "public: static void __cdecl Eigen::internal::general_matrix_matrix_product<__int64,double,1,0,double,0,0,0>::run(__int64,__int64,__int64,double const *,__int64,double const *,__int64,double *,__int64,double,class Eigen::internal::level3_blocking<double,double> &,struct Eigen::internal::GemmParallelInfo<__int64> *)" (?run@?$general_matrix_matrix_product@_JN$00$0A@N$0A@$0A@$0A@@internal@Eigen@@SAX_J00PEBN010PEAN0NAEAV?$level3_blocking@NN@23@PEAU?$GemmParallelInfo@_J@23@@Z)

Third,
It looks like Eigen matrix inversion doesnt utilize LAPACK functions.

For example I'm inverting next 3x3 double precision matrix A

7.1281933593749986499688020558096 4.8413655598958342807236476801336 -16.885416666666671403618238400668
4.8413655598958333925452279800083 3.2885666232638888217820749559905 -11.46875
-16.885416666666664298190880799666 -11.46875 40

Eigen inversion returns A_inv

3832.0164161361262813443318009377 1.0433858677166215814150343271175e-08 1617.6298464987853549246210604906
-1.0433858677166215814150343271175e-08 3832.0164161048246569407638162374 1098.7109567987947684741811826825
1617.6298464935682659415761008859 1098.7109568040118574572261422873 997.90488140786271742399549111724

A*A_inv =
0.99999998705607140436768531799316 9.19899321161210536956787109375e-08 -7.028575055301189422607421875e-09
2.947854227386415004730224609375e-08 1.0000000061354512581601738929749 -4.773028194904327392578125e-09
-1.028165570460259914398193359375e-07 -2.1790766812240897452433057487203e-07 0.99999996030601323582231998443604

My dumb lapack based implementation returns A_inv

3832.0166320546240967814810574055 1.1011253820762405399222080425344e-08 1617.6299376491328985139261931181
-8.14907252788543701171875e-10 3832.0166320172706946323160082102 1098.7110187108576155878836289048
1617.629937645741847518365830183 1098.7110187158500593795906752348 997.90493763799440785078331828117

A*A_inv =
0.99999999999636202119290828704834 0 0
-3.63797880709171295166015625e-12 1 0
1.4551915228366851806640625e-11 2.120868447542398278803446598495e-13 1

So lapack implementation is more preciese

Dumb LAPACK matrix inversion code:

Code: Select all

Eigen::Matrix3d inverse(const Eigen::Matrix3d &m) {
      std::vector<double> mci = { m(0, 0), m(0, 1), m(0, 2), m(1, 0), m(1, 1), m(1, 2), m(2, 0), m(2, 1), m(2, 2) };
      std::vector<int> ipiv(10);
      std::vector<double> work(9);
      int info, rows = 3, cols = 3, lwork = 9;
      int lda = 3;

      // LU factor
      dgetrf_(&rows, &cols, &mci[0], &lda, &ipiv[0], &info);

      // Invert
      dgetri_(&cols, &mci[0], &lda, &ipiv[0], &work[0], &lwork, &info);

      Eigen::Matrix3d minv;

      minv << mci[0], mci[1], mci[2], mci[3], mci[4], mci[5], mci[6], mci[7], mci[8];
      
      return minv;
   }


Hi,

four point 1), I'll quote LAPACKE.h file:

Code: Select all

* One can also redefine the types with his own types
* for example by including in the code definitions like
*
* #define lapack_complex_float std::complex<float>
* #define lapack_complex_double std::complex<double>
*
* or define these types in the command line:
*
* -Dlapack_complex_float="std::complex<float>"
* -Dlapack_complex_double="std::complex<double>"


We cannot really solve this issue on or side because the proper definition of these macro depend on how your lapack* libraries have been compiled.

So in your case you should define them before including any Eigen's header, instead of modifying lapacke.h.


For point 2), I cannot help much. What if you try yo call dgemm_ yourself similar to your inverse() example ? If that does not work neither, then you should check whether your BLAS.lib has been properly compiled with 'int' interface (and not 'long int') ? Also you could check whether BLAS.lib does contain a dgemm_ symbol by inspecting it: https://stackoverflow.com/questions/305 ... ibrary-lib


For point 3), it's simple: for compile-time 2x2, 3x3, and 4x4 matrices of float and double, Eigen uses special routines for much higher performance. You can enforce the default path using either a MatrixXd, or by enforcing the use of the LU factorization: Ainv = A.lu().inverse();