Intel® oneAPI Math Kernel Library
Ask questions and share information with other developers who use Intel® Math Kernel Library.
7093 Discussions

Problem with dstebz and dormtr (dense eigenvalue solver)

asd__asdqwe
Beginner
720 Views

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

0 Kudos
6 Replies
asd__asdqwe
Beginner
720 Views

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.) ?

0 Kudos
Ying_H_Intel
Employee
720 Views

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 ]

0 Kudos
asd__asdqwe
Beginner
720 Views

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.

0 Kudos
Aleksandr_Z_Intel
720 Views

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.

 

0 Kudos
Ying_H_Intel
Employee
720 Views

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  

0 Kudos
asd__asdqwe
Beginner
720 Views

Hello Ying,

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

0 Kudos
Reply