topic Re: Issue with workspace size generated by DGEEV in Intel® oneAPI Math Kernel Library

Issue with workspace size generated by DGEEV

david_engel — Wed, 22 Jan 2025 21:53:05 GMT

When using DGEEV in the below example, the first element of work (supposed to be the optimal workspace size for the calculation) is computed as 1. The call requires a minimum value of 8, so if we do a call to compute the workspace size and then do a call with that size it fails. Older versions of MKL as well as when calling SGEEV with the same inputs all return 68 as the optimal workspace size.

I first noticed this when moving from a very old version of MKL 2018, to the 2024 version, but it is still present in the latest version (INTEL_MKL_VERSION 20250001)

#include "mkl.h" typedef double dtype; int main() { const char* jobvl = "N"; const char* jobvr = "N"; const int n = 2; dtype a[] = {1,3,3,1}; const int lda = 2; dtype wr[] = {0,0}; dtype wi[] = {0,0}; const int ldvl = 2; const int ldvr = 2; dtype w = 8; dtype* work = &w; int lwork = 8; int info = -1; DGEEV( jobvl, jobvr, &n, a, &lda, wr, wi, NULL, &ldvl, NULL, &ldvr, work, &lwork, &info ); return work[0]; // returns 1 }

Re: Issue with workspace size generated by DGEEV

Sergey_K_Intel1 — Sat, 01 Feb 2025 01:53:12 GMT

There is a mistake in your code. Namely in order to get the optimal size of the work array, lwork (not info = -1 as you set) must be set to -1.

If we do this change (e.g. replace int lwork = 8; with int lwork = -1; ) , you code return 6 and 6 is the correct number according to the documentation (please take a look at the Lapack User's Guide https://www.netlib.org/lapack/lug/ )

In fact, the code you provided computes eigenvalue and it doesn't return the optimal size, and if we insert printing of the computed eigenvalue like this

printf( " first eigenvalue %15.8e +i %15.8e \n", wr[0], wi[0]);
printf( " second eigenvalue %15.8e +i %15.8e \n", wr[1], wi[1]);

we get

$ ./a.out

first eigenvalue 4.00000000e+00 + i 0.00000000e+00
second eigenvalue -2.00000000e+00 + i 0.00000000e+00

as it should be.

All the best

Sergey

Re: Issue with workspace size generated by DGEEV

david_engel — Mon, 03 Feb 2025 20:43:10 GMT

Hi Sergey,

I think there is some confusion here. If I supply -1 for lwork, then we return 6, which is a valid workspace size to be used in subsequent runs, but as per the ?geev documentation:

"If you are in doubt how much workspace to supply, use a generous value of lwork for the first run or set lwork = -1.

If you choose the first option and set any of admissible lwork sizes, which is no less than the minimal value described, the routine completes the task, though probably not so fast as with a recommended workspace, and provides the recommended workspace in the first element of the corresponding array work on exit. Use this value (work(1)) for subsequent runs."

So regardless of whether or not I set lwork to -1, this should complete the task, and then the optimal workspace size should be the first element of of the work array. The first element of the work array after the run is 1. Or if it means index 1, the result of work[1] is 3. For both of these values, if I use it as lwork in a subsequent run, the run fails.

Furthermore, if I do this exact same computation with SGEEV, it completes successfully and the first element of the work array is 68, which is the value it was in older versions of MKL.

Re: Issue with workspace size generated by DGEEV

Sergey_K_Intel1 — Mon, 03 Feb 2025 22:37:24 GMT

There are two different values of lwork used by all LAPACK routines and they behave differently depending on its value. First, all LAPACK routines always return the optimal size of the work array if lwork is set to -1. The optimal size depends on OneMKL or Reference LAPACK version and/or precision type and/or the type of OS and/or the number of threads and/or CPU. It’s not a constant and it may vary. This value is recommended for getting the best performance.

At the same time LAPACK routines are using the minimal required value of LWORK and online Reference LAPACK documentation always mentions the minimal value of LWORK. For example, Netlib Reference LAPACK for DGEEV says

“ LWORK is INTEGER

The dimension of the array WORK. LWORK >= max(1,3*N), and

if JOBVL = 'V' or JOBVR = 'V', LWORK >= 4*N. For good

performance, LWORK must generally be larger.”

Please also pay attention to the following sentence in OneMKL Developer Reference (https://www.intel.com/content/www/us/en/docs/onemkl/developer-reference-fortran/2025-0/geev.html

“Note that if you set lwork to less than the minimal required value and not -1, the routine exits immediately with an error and does not provide any information on the recommended workspace.“

So DGEEV as well as any another LAPACK must work if and only if LWORK is greater or equal to the minimal required value of LWORK (6 in your case because of jobr=’N’ jobvl=”N”). If LWORK is less than the minimal required value, you get an error message.

By the way you can use the LAPACKE interface for DGEEV (please see https://www.intel.com/content/www/us/en/docs/onemkl/developer-reference-c/2025-0/geev.html). In the case of using LAPACKE_dgeev you don’t need to take care of the LWORK size.

All the best

Sergey