Problem with dstebz and dormtr (dense eigenvalue solver)

[asd__asdqwe]

Hello,

I'm trying to solve a symmetric generalized eigenvalue problem with the BLAS Fortran bindings inside C and I'm facing two problem with the following piece of code:

1. if I want all eigenvalues (-DSEGFAULT), I get a segfault inside the call to dstebz,
2. I am not sure on how to call dormtr to get the correct eigenvectors of my original system.

This is a simple Neumann problem, so the first eigenvalue is 0, but the associated eigenvector should be constant and this is clearly not the case.

Thank you in advance for your help, I'm compiling with the following line: icpc gevp.cpp -L${MKLROOT}/lib/intel64 -lmkl_intel_lp64 -lmkl_core -lmkl_intel_thread -lmkl_rt -liomp5

[asd__asdqwe]

Concerning 1., it really is a MKL bug, since the same code when changing the include of "mkl.h" by "clapack.h" and compiled like this "icpc gevp.cpp -I/System/Library/Frameworks/vecLib.framework/Headers/ -framework vecLib -DSEGFAULT" does not segfault. Can you reproduce the problem on your side ? Could you help me for my real problem (2.) ?

[Ying_H_Intel]

Hello,

Right, we can produce the first problem with intel64. It seems to be an MKL bug, we are investigating it.  (And if i compile it with 32bit mkl,it seems works )

Regarding the function dormtr, there is Fotran sample code in  MKL install directory => example_core=>lapack=>dormtrx.f.

The call steps is like

Reduce A to tridiagonal form T = (Q**T)*A*Q
*
         CALL DSYTRD(UPLO,N,A,LDA,D,E,TAU,WORK,LWORK,INFO)
*
*        Calculate the two smallest eigenvalues of T (same as A)
*
         CALL DSTEBZ('I','B',N,VL,VU,1,4,ZERO,D,E,M,NSPLIT,W,IBLOCK,
     +               ISPLIT,WORK,IWORK,INFO)

         Calculate the eigenvectors of T, storing the result in Z
*
            CALL DSTEIN(N,D,E,M,W,IBLOCK,ISPLIT,Z,LDZ,WORK,IWORK,IFAILV,
     +                  INFO)
*

    CALL DORMTR('Left',UPLO,'No transpose',N,M,A,LDA,TAU,Z,
     +                     LDZ,WORK,LWORK,INFO)
*

I checked your C code, it looks the call should be no problem,  except  you alloced memoy with increasing way (is there any special reason for that?)  and a tiny problem on print part cout << " " << z[i+k*n];).   Mover,  as there is 0 as eigenvalue, the computing precision error was accumulated to last result.

Best Regards,

Ying

32bit result

Eigenvalues in ascending order:
4.56201e-16
2
3.33333
4

Eigenvector associated with eigenvalue 4.56201e-16:
V[0] = [ -0.67082 -0.447214 -0.447214 -0.387298 ]^T
A V[0] = 4.56201e-16 B * V[0] =>
[ -0.894427 ] = 4.56201e-16 * [ -1.72894 ]
[ -1.01426 ] = 4.56201e-16 * [ -0.894427 ]
[ -1.01426 ] = 4.56201e-16 * [ -0.894427 ]
[ -1.54919 ] = 4.56201e-16 * [ -0.774597 ]

Eigenvector associated with eigenvalue 2:
V[1] = [ -1.21024e-16 -0.707107 0.707107 5.55112e-17 ]^T
A V[1] = 2 B * V[1] =>
[ -2.22045e-16 ] = 2 * [ -1.86536e-16 ]
[ -2.82843 ] = 2 * [ -1.41421 ]
[ 2.82843 ] = 2 * [ 1.41421 ]
[ 2.22045e-16 ] = 2 * [ 1.11022e-16 ]

Eigenvector associated with eigenvalue 3.33333:
V[2] = [ 0.547723 -0.547723 -0.547723 0.316228 ]^T
A V[2] = 3.33333 B * V[2] =>
[ 4.38178 ] = 3.33333 * [ 1.41167 ]
[ -2.82335 ] = 3.33333 * [ -1.09545 ]
[ -2.82335 ] = 3.33333 * [ -1.09545 ]
[ 1.26491 ] = 3.33333 * [ 0.632456 ]

Eigenvector associated with eigenvalue 4:
V[3] = [ -0.5 0 -2.77556e-17 0.866025 ]^T
A V[3] = 4 B * V[3] =>
[ -2 ] = 4 * [ -0.133975 ]
[ -1.73205 ] = 4 * [ 0 ]
[ -1.73205 ] = 4 * [ -5.55112e-17 ]
[ 3.4641 ] = 4 * [ 1.73205 ]

[asd__asdqwe]

Hello Ying,

Thanks for your time. Ok, the bug is not critical for my application, but I guess it's better if you could fix it anyway. Concerning my example, I had forgotten to call dpotrs at the end to go back to the solution of the original generalized eigenvalue problem which was reduced by dsygst, but now the solution is correct, thanks ! :) The memory is allocated in an increasing way to avoid fragmentation.

[Aleksandr_Z_Intel]

Hi, qweasd q.

As a workaround you can try to pass something instead of NULL to dstebz:

    double vl, vu;

    dstebz_(&uplo, &order, &n, &vl, &vu, &il, &iu, &tol, d, e, &m, &nsplit, w, iblock, isplit, work, iwork, &info); 

But anyway the issue is on our side, we will let you know in this thread when it will be fixed.

[Ying_H_Intel]

Hello qweasd q, 

I'm glad to notify that MKL 11.1 update 4 and MKL 11.2 are available now. The problem should be fixed in the versions. You are welcomed to try them.

Best Regards,

Ying  

[asd__asdqwe]

Hello Ying,

Thank you for letting me know, I will try that ASAP.