Extended eigensolver problem

[Hainan_W_]

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;

}

[Gennady_F_Intel]

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.

[Hainan_W_]

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

[Irina_S_Intel]

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

[Hainan_W_]

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

[Irina_S_Intel]

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