topic feast works by the default in Intel® oneAPI Math Kernel Library

Extended eigensolver problem

Hainan_W_ — Fri, 12 Feb 2016 00:16:31 GMT

Hi all,

I am using MKL extended eigensolver (version 2015) for my eigenvalue problems. I tested it for a simple diagonal matrix {1, 2, 3}. However, I could not get the eigenvalues/eigenvectors using dfeast_scsrev. The returned info is "-4", seemingly to indicate that the matrix is not positive definite according to https://software.intel.com/en-us/node/470388#GUID-E1DB444D-B362-4DBF-A1DF-DA68F7FB7019.

I was not able to figure the issue with this simple problem. Can anyone help to take a look? The code was posted below.

Regards,

Hainan

#include <mkl.h>

bool test()
{
   int n =3;
   char uplo = 'U';
   int *iv, *jv;
   double *data;

   data = new double [3];
   iv = new int [4];
   jv = new int [3];

   data[0] = 1; data[1] = 2; data[2] = 3;
   iv[0] = 0; iv[1] = 1; iv[2] = 2; iv[3] = 3;//CSR format
   jv[0] = 0; jv[1] = 1; jv[2] = 2;

   //https://software.intel.com/en-us/node/470408
   //   dfeast_scsrev(const char* uplo , const MKL_INT* n , const double* sa , const MKL_INT* isa , const MKL_INT* jsa , MKL_INT* fpm , double* epsout , MKL_INT* loop ,
   //                   const double* emin , const double* emax , MKL_INT* m0 , double* e , double* x , MKL_INT* mode , double* res , MKL_INT* info);

   double epsout = 0.0;   //On output, contains the relative error on the trace: |tracei - tracei-1| /max(|emin|, |emax|)
   MKL_INT loop = 0;       //On output, contains the number of refinement loop executed. Ignored on input.
   double emin = -1000;       //The lower bounds of the interval to be searched for eigenvalues; emin < emax.
   double emax = 1000;       //The upper bounds of the interval to be searched for eigenvalues; emin < emax.
   MKL_INT m0 = 3;           //Total number of eigenvalues to be computed
   double *eigVal = new double[m0];
   double *eigVec = new double[m0*n];
   MKL_INT m = m0;           //The total number of eigenvalues found in the interval [emin, emax]: 0 <= m <= m0.
   double *res = new double[m0];
   MKL_INT info = 0;       //On output, if info=0, the execution is successful. If info != 0, see Output Eigensolver info Details.
   //https://software.intel.com/en-us/node/470388#GUID-E1DB444D-B362-4DBF-A1DF-DA68F7FB7019
   MKL_INT fpm[128];       //see https://software.intel.com/en-us/node/470386
   feastinit (fpm);       //Initialize Extended Eigensolver input parameters with default values. see https://software.intel.com/en-us/node/521737
   //fpm[0] = 1;               //Extended Eigensolver routines print runtime status to the screen.
   fpm[3] = 6;           //Error trace double precision stopping criteria

dfeast_scsrev(&uplo, &n, data, iv, jv, fpm, &epsout, &loop, &emin, &emax, &m0, eigVal, eigVec, &m, res, &info);

   delete data;
   delete iv;
   delete jv;
   delete eigVal;
   delete eigVec;
   delete res;

return true;

}

feast works by the default

Gennady_F_Intel — Fri, 12 Feb 2016 10:06:22 GMT

feast works by the default with one-based indexing, therefore, the following modification should help. please check it.

iv[0] = 1; iv[1] = 2; iv[2] = 3; iv[3] = 4;//CSR format
jv[0] = 1; jv[1] = 2; jv[2] = 3;

the error code == -4 in that case is the wrong output, actually it should be -104. this issue has been already addressed and will be fixed into the nearest 11.3 update 2 version of mkl.

Hi Gennady,

Hainan_W_ — Fri, 12 Feb 2016 20:34:58 GMT

Hi Gennady,

Thank you for looking into this problem. It indeed works after changing the index as you suggested. The resulting eigenvalues and eigenvectors are also close enough to theoretical values. Then, I also tested another symmetric matrix =

0 1 1

1 0 1

1 1 0

Its eigenvalues are -1, -1, 2 and the eigenvectors are {-1, 1, 0}, {-1, 0, 1}, {1, 1, 1}

(http://tutorial.math.lamar.edu/Classes/DE/LA_Eigen.aspx)

The input matrix in code is:

        data[0] = 1; data[1] = 1; data[2] = 1;
       iv[0] = 1; iv[1] = 3; iv[2] = 4; iv[3] = 4;//CSR format
       jv[0] = 2; jv[1] = 3; jv[2] = 3;

Now the feast solver gave correct eigenvalues but the corresponding eigenvectors look wired:

{-0.38573160677214430
-0.43035835023181784
0.81608995700396181}

{-0.71963726571353026
0.69387200336403154
0.025765262349498819}

{-0.57735026918962573
-0.57735026918962573
-0.57735026918962551}

Do you have any idea what the issue is?

Thanks!

Hainan

Hi Hainan,

Irina_S_Intel — Mon, 15 Feb 2016 05:08:19 GMT

Hi Hainan,

I understand your concern but there is no issue here. Eigenvectors are not unique.
Firstly, eigenvectors may vary by a factor as it happened in case of eigenvalue 2:

{-0.57735026918962573, -0.57735026918962573, -0.57735026918962551} = -0.577350269189625 * {1,1,1}

Secondly, if the eigenvalue has geometric multiplicity m then any m linearly independent vectors in an eigenspace correspondent to this eigenvalue can be chosen as eigenvectors. In your case, the geometric multiplicity of eigenvalue 1 is 2 and both sets { {-1, 1, 0}, {-1, 0, 1}} and{{-0.38573160677214430, -0.43035835023181784, 0.81608995700396181},{-0.71963726571353026, 0.69387200336403154,0.025765262349498819}} form a basis for the same eigenspace. It is easy to check that they are actually eigenvectors by verifying that Ax + 1*x = 0.

Best regards,

Irina

Hi Irina,

Hainan_W_ — Tue, 16 Feb 2016 06:41:00 GMT

Hi Irina,

Thank you for the reply! You are correct. I got it.

Now another question I have is: is there any way to solve for the largest couple of eigenvalues (by specifying an integer) for a matrix using mkl eigensolver (feast)? I seems the standard way is to specify a range of eigenvalues for the eigensolver (feast). But in my problem, I am most interested in eigenvalues with large magnitude (positive definite matrix). Moreover, it is hard to guess the range or magnitude of my eigenvalues. Therefore, specifying a range is inconvenient for me. Is there a better way to use eigensolver for my problem?

Thanks again!

Hainan

Irina S. (Intel) wrote:

Hi Hainan,

I understand your concern but there is no issue here. Eigenvectors are not unique.
Firstly, eigenvectors may vary by a factor as it happened in case of eigenvalue 2:

{-0.57735026918962573, -0.57735026918962573, -0.57735026918962551} = -0.577350269189625 * {1,1,1}

Secondly, if the eigenvalue has geometric multiplicity m then any m linearly independent vectors in an eigenspace correspondent to this eigenvalue can be chosen as eigenvectors. In your case, the geometric multiplicity of eigenvalue 1 is 2 and both sets { {-1, 1, 0}, {-1, 0, 1}} and{{-0.38573160677214430, -0.43035835023181784, 0.81608995700396181},{-0.71963726571353026, 0.69387200336403154,0.025765262349498819}} form a basis for the same eigenspace. It is easy to check that they are actually eigenvectors by verifying that Ax + 1*x = 0.

Best regards,

Irina

Hi Hainan,

Irina_S_Intel — Tue, 16 Feb 2016 10:18:57 GMT

Hi Hainan,

At the moment Extended Eigensolver doesn't support this operation. However, there are two ways how you can find k max or min eigenvalues using MKL functionality. First is to convert your matrix to dense storage format and use Lapack routines (https://software.intel.com/en-us/node/521045#E20540C0-34FD-4CD2-B803-B6096460B4F2). Second is to run mkl eigensolver feast on sufficiently large interval that for sure covers the largest eigenvalue. Algorithm will return info = 3 if this interval contains more than requested number of eigenvalues. Thus you can iteratively move left edge of an interval until it contains exactly k largest eigenvalues.

Best regards,

Irina