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

[Gael Guennebaud <gael.guennebaud@xxxxxxxxx>]

I see.

Problem fixed: https://bitbucket.org/eigen/eigen/commits/3a08a093adac (in 3.3 too).

gael

On Mon, Apr 9, 2018 at 9:52 AM, Peter <list@xxxxxxxxxxxxxxxxx> wrote:

Dear Gael,

to be sure I checked the eigen values I get after calling

eigensolverH.compute(  *itH->second )

and all samples I looked add where sorted.
Using nm the only LAPACK routine that is called is dsysev

000000000000752e t _GLOBAL__sub_I__ZN10TpOperator11DiagonalizeERN5Eigen6MatrixIdLin1ELi1ELi0ELin1ELi1EEERS_
                 U __gxx_personality_v0
                 U LAPACKE_dsyev

Am 07.04.2018 um 01:06 schrieb Gael Guennebaud:
> It's not clear to me what you call "mess". Do you get completely different eigenvalues for the same input matrix? Or do your iterations slightly diverge before exploding?

I get a few correct than I get complete nonsense. I think I found the problem, see below.

Looking at the description it says:

< https://software.intel.com/en-us/node/521120 >

    On exit, if jobz = 'V', then if info = 0, array a contains the orthonormal eigenvectors of the matrix A.

    If jobz = 'N', then on exit the lower triangle

    (if uplo = 'L') or the upper triangle (if uplo = 'U') of A, including the diagonal, is overwritten.

w: If info = 0, contains the eigenvalues of the matrix A in ascending order.

So, indeed the eigen values should be sorted.
I tried to diagonalize a copy of the matrix, to check whether something gets overwritten,

                        TpMatrix HHH = *itH->second;
                        eigensolverH.compute(  HHH );

but that doesn't change something.

However changing

U( itH->first.first , itH->first.second, itH->second ->rows(), itH->second ->cols()   ) =  eigensolverH.eigenvectors();

to

U( itH->first.first , itH->first.second, itH->second ->rows(), itH->second ->cols()   ) =  eigensolverH.eigenvectors().adjoint();

the results look fine. In return I suspect that the MKL bindings provide a transposed transformation matrix.
The matrices are defines as

typedef Eigen::Matrix< double,   Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> TpMatrix;


Switching to ColMajor, the MKL and non-MKL version both need the first, non-adjoint, version.

Best regards,
Peter