topic Hi Areg, in Intel® oneAPI Math Kernel Library

Feast can't find eigenvalues in the specified interval

Areg_G_ — Tue, 19 Jul 2016 22:59:50 GMT

Hi,

I am trying to solve eigenvalue problem using zfeast_hcsrev routine, here is the excerpt from my code

MKL_INT fpm[128];
feastinit(fpm);
fpm[0] = 1;
fpm[1] = 20;
char uplo = 'U';
MKL_INT n = 85;
double emin = 11.4;
double emax = 12;
double epsout;
MKL_INT loop;
MKL_INT m0 = 20;
MKL_INT m;
double en[m0];
std::complex<double>* xeig = new complex<double>[nimbsize * m0];
double res[m0];
MKL_INT info;
zfeast_hcsrev(&uplo, &n, (const MKL_Complex16*) baseHamFull.mat, (const long long int*) baseHamFull.ia,
(const long long int*) baseHamFull.ja, fpm, &epsout, &loop, &emin, &emax,
&m0, en, (MKL_Complex16*) xeig, &m, res, &info);

And here is what I get:

Intel MKL Extended Eigensolvers: complex double precision driver
Intel MKL Extended Eigensolvers: List of input parameters fpm(1:64)-- if different from default
Intel MKL Extended Eigensolvers: fpm(1)=1
Intel MKL Extended Eigensolvers: fpm(2)=20
Search interval [1.140000000000000e+01;1.200000000000000e+01]
Intel MKL Extended Eigensolvers: Size subspace 20
#Loop | #Eig | Trace | Error-Trace | Max-Residual
Intel MKL Extended Eigensolvers WARNING: No eigenvalue has been found in the proposed search interval.
==>INFO code =: 1
Intel MKL Extended Eigensolvers WARNING: No eigenvalue has been found in the proposed search interval.
==>INFO code =: 1
1 0

You can see in the files the matrix which I have supplied in CSR format. I know that there are eigenvalues in this range, by solving with different library. Some of the eigenvalues are

11.8341 11.839 11.877 11.8801 11.8922 11.9136 11.9127 11.9215 11.929 11.9321 11.9489 11.9538 11.9518 11.9586 11.9678 11.9763 11.9837 11.9935 11.9974 12.0009

Can you say what I am doing wrong? Also in my problems usually I need lowest eigenvalues, but don't usually know the range where to look for them. Is there a way to implement this in Feast or not?

Thank you for your help.

Areg Ghazaryan

Hi Areg,

Irina_S_Intel — Tue, 26 Jul 2016 11:33:39 GMT

Hi Areg,

You can do a simple check to prove that the eigenvalues are indeed the ones that you state using Gershgorin circles. Gershgorin theorem states that all eigenvalues lie within so-called Gershgorin circles, thus we can calculate the interval that must contain all eigenvalues as the most left and most right intersection of these circles with the real line.

The following C code will do this, [min, max] is an interval that contains all eigenvalues for sure:

void gershgorin_crcl(int n, int *ia, int *ja, MKL_Complex16 *a, double *min, double* max)
{
    int row, col, i, j;
    double *rad = (double *)malloc( n * sizeof(double) );
    int start_row;....
    int end_row;
    double left_edge, right_edge;
    *min = 100000000.0;
    *max = -100000000.0;

    for( i = 0; i<n; i++ )
    {
<------>rad = 0.0;
    }

    for( row = 0; row<n; row++ )
    {
<------>start_row = ia[row] - 1;
<------>end_row   = ia[row+1] - 1;
<------>for (j = start_row + 1; j<end_row; j++)
<------>{
<------>    col = ja-1;
<------>    double abs_a = sqrt(a.real*a.real+a.imag*a.imag);
<------>    rad[row] += abs_a;
<------>    rad[col] += abs_a;
<------>}
    }
    for ( row = 0; row<n; row++ )
    {
<------>start_row = ia[row] - 1;
<------>left_edge = a[start_row].real - rad[row];
<------>right_edge = a[start_row].real + rad[row];
<------>if ( *min > left_edge ) *min = left_edge;
<------>if ( *max < right_edge )*max = right_edge;<---->
    }
    free(rad);
}

I ran it on your matrix and the answer I got is the following:

min = 497.979941 max = 580.954978

Now I run zfeast_hcsrev on [min, max] and get all 85 eigenvalues.

E[0]=525.897983
E[1]=526.117500
E[2]=527.806863
E[3]=527.943623
E[4]=528.479884
E[5]=529.392685
E[6]=529.433090
E[7]=529.781798
...

E[82]=550.344795
E[83]=553.120385
E[84]=556.137079

Please make sure that your matrix is correct.

Best regards,

Irina

Hi Irina,

Areg_G_ — Tue, 26 Jul 2016 20:52:13 GMT

Hi Irina,

Thank you for your help. I had indeed an error, in the sense that in different package I had multiplication factor, so the range which I was specifying for Feast was wrong. Once I have corrected that I got the right eigenvalues. I have two separate questions:

1. I noticed that if I increase emax, so that m gets closer to m0, the found eigenvalues become wrong. I know that this was written in the documentation, so this is natural. So do you have some kind of suggestion how to specify emin and emax. The Gershgorin method, I think is not helpful for me, because the matrix which I have provided is just an example. Usually my matrices have a size from 500000 to 1000000 and I only need smallest 30-40 eigenvalues, therefore specifying emax is still unclear for me, because the eigenvalues are closely spaced and if I specify slightly bigger range I get wrong result. Of course I can write iterative code, so that if m0 < 1.5 m to decrease emax, but that will make the calculation very time consuming, because it should get the eigenvalues several times.

2. My second question is related to the compiler. I compiled first locally at my computer with MKL 11.3.2 and gcc 4.8 and I got the results. Later I tried to compile on shared memory machine with mkl/2015.1.135 and intel compiler 15.0.1.133, but I get segmentation fault. Here is my compilation command:

icpc job101.cpp basis.cpp opermask.cpp -std=c++11 -Wall -g -DMKL_ILP64 -qopenmp -mkl=parallel -o job101 -L/share/apps/gsl/1.16/lib/ -lpthread -lm -ldl -lgsl -lgslcblas

What I am doing wrong here.

Thank you.

Areg

Hi Areg,

Irina_S_Intel — Wed, 27 Jul 2016 06:08:00 GMT

Hi Areg,

1. Wrong eigenvalues is not expected behavior of Extended Eigensolver, can you please provide a reproducer of the issue? We indeed recommend to set m0 a 1.5 times larger but that is to make sure that all eigenvalues within the search interval are found. Eigenvalues that lie very close to the edge of an interval can probably be lost but all the eigenvalues found will be correct. For that reason you can always make your interval a bit larger. There is also a way how you can calculate number of eigenvalues within your interval exactly using Sylvester's law of inertias:

run phase 11, phase 22 of Intel MKL PARDISO with matrix A-emaxI (I - identity matrix) with iparm[22]=1. After phase 22 iparm[22] will contain number of negative eigenvalues of A-emaxI or (what's the same) number of eigenvalues of A that are smaller than emax. This will give you an exact number of eigenvalues within [emin, emax]. For more info please see Pardiso documentation page https://software.intel.com/en-us/node/521678.

2. Looks fine, can you please send me source files so I could reproduce it with your options.

Best regards,

Irina

Hi Irina,

Areg_G_ — Thu, 28 Jul 2016 00:14:32 GMT

Hi Irina,

1. Please see the source code attached. Change emax on line 93 in job101.cpp 780, 785, 790 and you will see that each time it will give different eigenvalues. Please check that out, and if it is normal then I will use Sylvester's law of inertias, which you have suggested.

2. I compile the code with gcc

g++ job101.cpp basis.cpp opermask.cpp -std=c++11 -Wall -O2 -DMKL_ILP64 -fopenmp -m64 -o job101 -L/share/apps/gsl/1.16_gcc/lib -lgsl -lgslcblas -Wl,--no-as-needed -L${MKLROOT}/lib/intel64 -lmkl_intel_ilp64 -lmkl_core -lmkl_gnu_thread -lpthread -lm

and with intel compiler

icpc job101.cpp basis.cpp opermask.cpp -std=c++11 -Wall -g -DMKL_ILP64 -qopenmp -mkl=parallel -o job101 -L/share/apps/gsl/1.16/lib/ -lpthread -lm -ldl -lgsl -lgslcblas

As I said before it runs normal when compiled with gcc, but gives segmentation fault when compiled with intel compiler. I presume, that this happens, because the old intel compiler (mkl version 2015.1.135 and intel compiler 15.0.1.133) does not fully use 64 bit integers in the internals of feast.

3. I have another issue related to memory. The memory consumption of feast is huge, for the matrix which uses 600 MB of memory, feast uses 50Gb. Is this normal or something should be changed in the code?

Thanks,

Areg

Hi Areg,

Irina_S_Intel — Mon, 01 Aug 2016 07:34:00 GMT

Hi Areg,

1. The thing here is that you forgot to check info after feast run. In case of search interval [773, 785] zfeast_hcsrev returned info=3 that corresponds to m0 too small. In fact, there are 144 eigenvalues in that interval.

2. Unfortunately, this issue I can't reproduce on my side. Both intel and gcc compiler built codes work similar with no crashes. Old Intel compiler as well as feast support 64 bit integers.

3. That is not expected behavior. Can you give me a reproducer?

By the way, may I ask you from what kind of application is the eigenvalue problem that you tackle?

Best regards,

Irina

Hi

kumar__aman — Fri, 13 Jul 2018 16:49:00 GMT

Is there any reproducer? You

Gennady_F_Intel — Mon, 12 Nov 2018 03:45:47 GMT

Is there any reproducer? You may also try the latest mkl 2019 u1 we released the last week and check if the problem exists with this update.