Re: [eigen] problem with MKL bindings for dense Eigen

[Gael Guennebaud <gael.guennebaud@xxxxxxxxx>]

oh, one more thing, I'm not sure that the eigenvalues are ordered the same way when calling MKL, just in case you relied on that particular ordering.

gael

On Fri, Apr 6, 2018 at 10:49 PM, Gael Guennebaud <gael.guennebaud@xxxxxxxxx> wrote:

Just a guess, but maybe it's just that MKL gives you different values for the lowest eigenvalues (e..g., 1e-16 vs 1e-30 vs exact 0 while the largest is 1) but your code is lacking some thresholds to properly deal with them (e.g., to ignore them). So tracking the largest and the smallest eigenvalues might a be a good start.

gael

On Fri, Apr 6, 2018 at 9:14 PM, Peter <list@xxxxxxxxxxxxxxxxx> wrote:

Dear All,

the discussion on MKL and sparse Eigen routines reminded me,
that I used to have problems with the MKL bindings, so just tried it again.

I'm currently investigating a system, where a phase transition
including scaling laws is governed by the precision of the underlying
arithmetic. From a numerical point of view it is a nice testbed for the
precision of the libraries. In case you are curious, there's a preprint
on the arxiv server:  < https://arxiv.org/abs/1803.08480 > .

Taking the code for double precision and compiling with  -DEIGEN_USE_BLAS
I get similar results, as without MKL. Albeit I see small deviations already
in the stable regime. IN this regime, the MKL version appears to be more accurate
compared to plain Eigen. However, in the precision sensitive region, MKL is less accurate.


My real problem is, that compiling with  -DEIGEN_USE_BLAS -DEIGEN_USE_LAPACKE_STRICT
gives me a few iteration steps that are still correct and then complete nonsense.

The only LAPACK routine I'm aware of that could be called from the bindings should
be a diagonalization (yeah, the code is written by me, but some parts are a few years old):

Eigen::SelfAdjointEigenSolver<TpMatrix> eigensolverH;
.....
eigensolverH.compute(  *itH->second );  ///< within a loop over matrices

Before I spend time in hunting the problem, I'd like to ask whether there's
more to do than the above defines and the correct linking options.
Note, I'm using a 2018 MKL and compiled with g++ 7.2 on Linux 4.15.13

CFLAGS = -std=c++11 -O3 -march=native -DNDEBUG

linking with MKL is done via
MKLLIBS = -L $(HOME)/intel/MKL2018/mkl/lib/intel64 -lmkl_intel_lp64 -lmkl_intel_thread -lmkl_core -liomp5 -lpthread -lm -ldl
or
MKLLIBS = -L $(HOME)/intel/MKL2018/mkl/lib/intel64 -lmkl_intel_lp64 -lmkl_sequential -lmkl_core -lpthread -lm -ldl

both versions show the same nonsense.

Thanks in advance,
best regards,
Peter