MKL Eigensolvers for Fortran code

[Ioannis_K_]

Hello,

I am interested in checking the available MKL eigensolvers in a fortran program I have, both for linear and generalized eigenvalue problems.

I want to be able to calculate the 'k' lowest eigenvalues of the problem at hand.

I would be grateful if it was possible for someone to provide me with a summary of what functions (feast etc.) are available on mkl for this type of problems. I have already found the feast mkl documentation, but I understand there is also a Krylov-Schur method available in mkl.

Furthermore, the feast documentation indicates that it can be used to find the eigenvalues within an interval, but not how to find the 'k' lowest eigenvalues.

Thank you in advance.

[VidyalathaB_Intel]

Hi Ioannis,

Thanks for reaching out to us.

Could you please let us know if ?spevx routine (which computes selected eigenvalues) matches your use case in finding the 'k' lowest eigenvalues ?

Please refer to the below link for more details regarding the same

https://www.intel.com/content/www/us/en/develop/documentation/onemkl-developer-reference-fortran/top/lapack-routines/lapack-least-squares-and-eigenvalue-problem/lapack-least-squares-eigenvalue-problem-driver/symmetric-eigenvalue-problems-lapack-driver/spevx.html

Also here is the list of LAPACK driver routines used for solving symmetric eigenvalue problems.

Please do let us know if this is not what you are looking for.

Regards,

Vidya.

[Ioannis_K_]

Hi Vidya,

Thank you for your reply and for all this information.

I had one follow-up question: is there a function specifically tailored to sparse matrices? I have very large, symmetric matrices in CSR format, so I would prefer to use something like feast (which is specifically tailored to CSR). If not, is there an mkl function for reducing a sparse matrix to, e.g., tridiagonal form, so that I can subsequently use one of the LAPACK functions for finding eigenvalues/eigenvectors?

Again, thank you for your help.

Yannis

[VidyalathaB_Intel]

Hi Ioannis,

Thanks for getting back to us.

Could you please let us know if you have gone through the function mkl_sparse_?_gv https://www.intel.com/content/www/us/en/develop/documentation/onemkl-developer-reference-fortran/top/extended-eigensolver-routines/extended-eigensolver-interfaces/ext-eigensolve-ifaces-find-large-small-eigenvalues/mkl-sparse-gv.html ?

Meanwhile, you can refer to the below links for some functions which might help you

https://www.intel.com/content/www/us/en/develop/documentation/onemkl-developer-reference-c/top/blas-and-sparse-blas-routines/sparse-blas-level-2-and-level-3-routines/sparse-blas-level-2-and-level-3-routines-1/mkl-csrdia.html

https://www.intel.com/content/www/us/en/develop/documentation/onemkl-developer-reference-c/top/lapack-routines/lapack-least-squares-and-eigenvalue-problem/lapack-least-square-eigenvalue-problem-computation/symmetric-eigenvalue-problems-lapack-computation/sptrd.html

Regards,

Vidya.

[VidyalathaB_Intel]

Hi,

As we haven't heard back from you, could you please provide us with an update regarding the issue?

Regards,

Vidya.

[VidyalathaB_Intel]

Hi,

We assume that your issue is resolved and closing this thread. Please post a new question if you need any additional assistance from Intel as this thread will no longer be monitored.

Regards,

Vidya.

[Ioannis_K_]

Hello Vidya,

First, I wanted to thank you for your previous suggestions in using mkl_sparse_d_gv() to calculate the lowest eigenvalues in a generalized eigenvalue problem.

I finally found time to work on this. However, when I call the function in my program, using a statement:

stat = mkl_sparse_d_gv(...) (I include the various arguments etc.)

In the above, stat is an integer variable whose value informs me whether the function successfully conducted the computation.

When I run it, the program does not compute any eigenvalues/vectors, and simply returns a value of stat = 5. 

I wanted to ask: is there a resource for me to check what is the meaning of various values of the integer stat? I would imagine that stat = 0 would imply a normal completion of the computation, but what is the meaning of stat = 5 or stat = 6 etc.?

Again, thank you very much for your help.

Yannis

[mecej4]

The symbolic constants returned by MKL routines are defined in the relevant include files in the MKL include directory. You can see the following definitions in mkl_spblas.h:

SPARSE_STATUS_SUCCESS           = 0,    /* the operation was successful */
SPARSE_STATUS_NOT_INITIALIZED   = 1,    /* empty handle or matrix arrays */
SPARSE_STATUS_ALLOC_FAILED      = 2,    /* internal error: memory allocation failed */
SPARSE_STATUS_INVALID_VALUE     = 3,    /* invalid input value */
SPARSE_STATUS_EXECUTION_FAILED  = 4,    /* e.g. 0-diagonal element for triangular solver, etc. */
SPARSE_STATUS_INTERNAL_ERROR    = 5,    /* internal error */
SPARSE_STATUS_NOT_SUPPORTED     = 6     /* e.g. operation for double precision doesn't support other types */

If your example code produced a status of "Internal Error", it will be necessary to post your source code to enable the MKL support staff to respond.

[Ioannis_K_]

Thank you very much for the reply. I apologize, I had not noticed that the documentation for mkl_sparse_d_gv(...) was actually explaining the error codes (I was looking for values of the integer stat, while the documentation was simply listing the "enumerated return values".

I was able to find the reason for my problem. I had to call the mkl_sparse_ee_init () function to initialize the integer parameter vector pm(), before passing it in the mkl_sparse_d_gv(...) function.

I had two follow-up questions: 

First, I want to solve a generalized eigenvalue problem: [A]{X}-lambda[B]{X} = {0}. The matrix [B] in my case is DIAGONAL (but non-unit). 

Right now, when I establish the MATRIX_DESCR structure for matrix [B], which I name Bdesc, I use the following:

Bdesc%type = SPARSE_MATRIX_TYPE_SYMMETRIC
Bdesc%mode = SPARSE_FILL_MODE_UPPER
Bdesc%diag = SPARSE_DIAG_NON_UNIT

Obviously, I also have to define the other integer arrays corresponding to the csr representation of [B]. 

Is there another way define the parameters of the MATRIX_DESCR structure, so that the code knows that [B] is a diagonal matrix?

The second question is: I am now able to run my code, but the performance of the algorithm is not as good as I would expect. I am running it for a large problem with a matrix dimension of about 55,000, and I try to find the lowest 20 eigenvalues. I build my code ensuring that I rely on the parallel version of mkl.

1) First of all, the Krylov-Schur algorithm (pm(3) = 1 - it is apparently also used when I set pm(3) = 0) is very slow in debug configuration. It takes about 12 minutes to complete (I have verified that it gives the correct lowest mode vectors and eigenvalues). I tried to relax the convergence tolerance (by changing the pm(2) parameter value), but this does not appear to affect the performance.

Even worse, when I run my code with the same exact problem in Release configuration (and optimization to maximize speed), the algorithm seems to be taking forever and has not even completed after 35 minutes! I had to stop the program, as I was not sure if and when it would complete. This alone is extremely strange, as I have never encountered situations where the optimized, release configuration of a program becomes so much slower than the debug configuration.

2) When I try using the subspace iteration algorithm (pm(3)=2), it never works; instead, it aborts with an internal error (stat = 5). This happens no matter how many eigenvectors I request (I tried 20, 50, and 500). 

Any advice on what might be the reason for the problem would be most appreciated. Is there something I could change in the pm() vector to speed up the Krylov-Schur algorithm or help the subspace iteration algorithm actually work? I looked at the pertinent documentation, but did not see many options for pm() values.