Experimenting with the symetric general eigen value solver dfeast_scsrgv.
Running it with a test matrix, I find that it overwrites the values for e and x on the stack?
dfeast_scsrgv (&uplo, &n, a, rows, columns,
b, rows, columns,
fpm, &epsout, &loop, &emin, &emax, &m0, eigs, eigv, &m, &res, &info);
exits with info=0, m=1, suggesting one eigenvalue was found but eigs and eigv have been changed and no longer point to the original arrays.
inputs are declared as
double* a = (double*) mkl_malloc(nnz*sizeof(double), MKL_ALIGN);
double* b = (double*) mkl_malloc(nnz*sizeof(double), MKL_ALIGN);
int* columns = (int*) mkl_malloc(nnz*sizeof(int), MKL_ALIGN);
int* rows = (int*) mkl_malloc((n+1)*sizeof(int), MKL_ALIGN);
char uplo = 'U'; //upper triangular matrices
MKL_INT m0 = 60;
double emin = 2;
double emax = 150;
double epsout = 0, res = 0;
double *eigs = (double*) mkl_malloc(m0 * sizeof(double), MKL_ALIGN);
double *eigv = (double*) mkl_malloc(m0 * n * sizeof(double), MKL_ALIGN);
MKL_INT loop=0, m=0, info = 0;
Setting a data breakpoint suggests it happens at mkl_core.dll!000007fee5371abf()
I thought I had the wrong size for one of the arguments, but they look ok, and there doesn't appear to be a variable adjacent on the stack.
Test project is attached, VS2010 with MKL 11.0.3
Anyone seen anything similar?
I find that it overwrites the values for e and x on the stack
I could not duplicate this behavior with the current version of the compiler (13.1.171) and MKL 11.0.3. I added statements to print the pointer values of eigs and eigv before and after the call to scsrgv, and set fpm=1 after the call to feastinit to enable more verbose output. I found that the pointer values were preserved; it is possible that the pointer values were saved, changed and restored within MKL, and a data breakpoint would catch that; however, I don't see why the library should not temporarily alter arguments passed to it as long as values that are expected to be preserved are restored before returning to the point of the call.
Note that a single eigenvalue is found, so your specifying m0=60 causes the array eigv to consume much more memory than is needed.
Thanks for taking a look.
Found the problem: I had res being declared as a double, instead of double[m0], so the solver writing to res was corrupting the stack for i>0. Missed it in the documentation initially, as everything gets passed in as a pointer.
Alan Richie wrote:
I had res being declared as a double, instead of double[m0]
Yes, that is an error which should be corrected. However, with the specific matrix in your example, there is only one eigenvalue in the range provided, so this problem would probably remain hidden since (i) passing the address of a scalar variable and (ii)the array name of an array with one element are.equivalent.