- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

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:

- if I want all eigenvalues (-DSEGFAULT), I get a segfault inside the call to dstebz,
- 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

Link Copied

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

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

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

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 ]

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

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.

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

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.

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

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

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

Hello Ying,

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

- Subscribe to RSS Feed
- Mark Topic as New
- Mark Topic as Read
- Float this Topic for Current User
- Bookmark
- Subscribe
- Printer Friendly Page