Intel MKL Feast for shift-invert diagonalization

[DebasishBanerjee]

Hi,

 I have a standalone code that uses FEAST to find all the eigenvalues of a matrix (in CSR format).

I am attaching the code with this email. I want to modify the code to find a fixed number (say 20) eigenvalues around a certain eigenvalue (E0) which I will specify in the code. Could you please suggest how this modification should look like?

 The documentation was not clear on this (namely on what parameters needed to be modified); and moreover I was unable to find any example which implements a similar problem.

I am attaching both the code and the input matrix (please untar + unzip it) in case anyone would like to try.

Best regards,
Debasish

=============================================================

#include <iostream>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "mkl.h"
#include "mkl_solvers_ee.h"
#include "mkl_types.h"
#include <fstream>
#include <vector>

int main()
{
/* Matrix A of size N in CSR format. Size N and 3 arrays are used to store matrix in CSR format */
int alignment=64;
unsigned int k;
/* read matrix A from a binary file */
std::ifstream inFile("sparse_mat.bin", std::ios::in | std::ios::binary);
int nrows, ncols, nval;
inFile.read(reinterpret_cast<char*>(&nrows), sizeof(int));
inFile.read(reinterpret_cast<char*>(&ncols), sizeof(int));
inFile.read(reinterpret_cast<char*>(&nval), sizeof(int));

std::cout<<"No of rows = "<<nrows<<std::endl;
std::cout<<"No of cols = "<<ncols<<std::endl;
std::cout<<"No of Hamiltonian entries = "<<nval<<std::endl;

std::vector<MKL_INT> rows(nrows);
std::vector<MKL_INT> cols(ncols);
std::vector<double> val(nval);
inFile.read(reinterpret_cast<char*>(rows.data()), nrows * sizeof(MKL_INT));
inFile.read(reinterpret_cast<char*>(cols.data()), ncols * sizeof(MKL_INT));
inFile.read(reinterpret_cast<char*>(val.data()), nval * sizeof(double));
inFile.close();

printf("size of rows = %d \n",nrows);
for(k=0;k<nrows;k++) printf("row[%d] = %d \n",k,rows[k]);
printf("size of cols = %d \n",ncols);
for(k=0;k<ncols;k++) printf("col[%d] = %d \n",k,cols[k]);
printf("elements = %d \n",nval);
for(k=0;k<nval;k++) printf("val[%d] = %f \n",k,val[k]);

char UPLO = 'F'; /* Type of matrix: (F means full matrix, L/U - lower/upper triangular part of matrix) */
const MKL_INT N = nrows-1;
/* Declaration of FEAST variables */
MKL_INT fpm[128]; /* Array to pass parameters to Intel MKL Extended Eigensolvers */
MKL_INT* r=&rows[0];
MKL_INT* c=&cols[0];
double* v=&val[0];
double Emin, Emax; /* Lower/upper bound of search interval [Emin,Emax] */
double epsout; /* Relative error on the trace */
MKL_INT loop; /* Number of refinement loop */
MKL_INT L = N; /* Should not be more than the dimension of the matrix */
MKL_INT M0; /* Initial guess for subspace dimension to be used */
MKL_INT M; /* Total number of eigenvalues found in the interval */
double *E; /* Eigenvalues */
double *X; /* Eigenvectors */
double *res; /* Residual */
E = (double*)(mkl_calloc(N,sizeof(double),alignment));
X = (double*)(mkl_calloc(N*N,sizeof(double),alignment));
res = (double*)(mkl_calloc(N,sizeof(double),alignment));

/* Declaration of local variables */
MKL_INT info; /* Errors */
double one = 1.0; /* alpha parameter for GEMM */
double zero = 0.0; /* beta parameter for GEMM */
MKL_INT i, j;
double trace, smax, eigabs;

printf("\n FEAST DFEAST_SCSREV AND DFEAST_SCSRGV EXAMPLE\n");
/* Initialize matrix X */
for (i=0; i<N*N; i++) X[i] = zero;

printf("Sparse matrix size %i\n", (int)N);

/* Search interval [Emin,Emax] */
Emin = -50.0;
Emax = 50.0;
printf("Search interval [ %.15e, %.15e ] \n", Emin, Emax);

M0 = L;
M = L;
loop = 0;
info = 0;
epsout = 0.0;

/* Step 1. Call FEASTINIT to define the default values for the input FEAST parameters */
feastinit(
fpm /* OUT: Array is used to pass parameters to Intel MKL Extended Eigensolvers */
);

fpm[0] = 1; /* Extended Eigensolver routines print runtime status to the screen. */

/* Step 2. Solve the standard Ax = ex eigenvalue problem. */
printf("Testing dfeast_scsrev routine:\n");
dfeast_scsrev(
&UPLO, /* IN: UPLO = 'F', stores the full matrix */
&N, /* IN: Size of the problem */
v, /* IN: CSR matrix A, values of non-zero elements */
r, /* IN: CSR matrix A, index of the first non-zero element in row */
c, /* IN: CSR matrix A, columns indices for each non-zero element */
fpm, /* IN/OUT: Array is used to pass parameters to Intel MKL Extended Eigensolvers */
&epsout, /* OUT: Relative error of on the trace */
&loop, /* OUT: Contains the number of refinement loop executed */
&Emin, /* IN: Lower bound of search interval */
&Emax, /* IN: Upper bound of search interval */
&M0, /* IN: The initial guess for subspace dimension to be used. */
E, /* OUT: The first M entries of Eigenvalues */
X, /* IN/OUT: The first M entries of Eigenvectors */
&M, /* OUT: The total number of eigenvalues found in the interval */
res, /* OUT: The first M components contain the relative residual vector */
&info /* OUT: Error code */
);
printf("FEAST OUTPUT INFO %d \n",info);
if ( info != 0 )
{
printf("Routine dfeast_scsrev returns code of ERROR: %i", (int)info);
return 1;
}
printf("Number of eigenvalues found %d \n", M);
for (i=0; i<M; i++){
printf("%.15e \n", E[i]);
}

return 0;
}

[Aleksandra_K]

Hi,

The function you are using computes all the eigenvalues within a given search interval (?feast_scsrev/?feast_hcsrev). There is another function mkl_sparse_?_ev which computes the largest/smallest eigenvalues of sparse matrices and you can specify the number of eigenvalues to be found as a parameter k0.  

Regards,

Alex

[DebasishBanerjee]

Hi,

 Could you please provide me with a sample code to use the function mkl_sparse_?_ev ?
 Moreover, I want to access a *few* eigenvectors in the middle of the spectrum but not at the edges, and would ideally like to apply the Shift-Invert method.

 I am experiencing quite a few problems with MKL Feast, with the sparse matrix the diagonalization often crashes (ie, a segmentation fault) or sometimes gets stuck (does not return from the function).

 Best,

 Debasish

[Aleksandra_K]

Hi, 

For the code examples please check the directory /opt/intel/oneapi/mkl/latest/share/doc/mkl/examples/c/sparse_eigsolvers/source. There are a few examples using the Feast algorithm and the file mkl_sparse_s_ev.c demonstrates the use of mkl_sparse_?_ev function. 

[Aleksandra_K]

Did you have a chance to try the suggested examples?

[DebasishBanerjee]

Hi,

 I have been trying the FEAST algorithm, but there seems to be a problem while passing the CSR matrix to the solver.

In several cases, the FEAST algorithm either crashes, or keeps running in an infinite loop without giving any error message. If it is helpful, I can share some detailed examples.

Best,

Debasish

[Aleksandra_K]

Please share some code example.

[DebasishBanerjee]

Hi Aleksandra,

 I am attaching the standalone cpp code as well as the binary file which I was trying it on.

 The binary file is attached in a gzipped file and uploaded here.

 I can confirm that:

 1. The supplied matrix is real symmetric and stored in a CSR format; it's size is 1107 x 1107 and has

     9651 non-zero entries. When I try to run the code, this is the output I get:

   ========================================================

No of rows = 1108
No of cols = 9651
No of Hamiltonian entries = 9651
size of rows = 1108
size of cols = 9651
elements = 9651

FEAST DFEAST_SCSREV AND DFEAST_SCSRGV EXAMPLE
Sparse matrix size 1107
Search interval [ -1.000000000000000e+01, 1.000000000000000e+01 ]
Testing dfeast_scsrev routine:

 ====================================================

   It clearly shows that binary matrix is read correctly, amd maybe there is some mistake in invoking
  the FEAST solver. Could you please help me with that? 

 2. A dense eigensolver has no problem to diagonalize this (such as the intel MKL lapack dsyev solvers).

 3. I really only a few eigenvalues, *but from the bulk of the spectrum*, using the shift-invert method.
     Consequently, I think the MKL feast eigenvolvers are ideal. The solver you suggested mkl_sparse_?_ev
     seem to only give eigenvalues/vectors from the edge of the spectrum, as stated in the documentation:

      https://www.intel.com/content/www/us/en/docs/onemkl/developer-reference-c/2025-1/mkl-sparse-ev.html

    If I am mistaken, please give me an example of how I can use the suggested solver to find eigenvalues in
    a band of higher energies.

  4. This is a test example. The actual problem will be a
        Matrix shape: (616227, 616227)
        Matrix nnz: 8644667
      Further I can confirm that using scipy's eigsh, I am easily able to get ~40-50 eigenvalues at the boundary,
      but eigsh (which wraps around ARPACK) has problems to find interior eigenvalues. This is however, known
      that ARPACK encounters problems. I am hoping that MKL feast would be able to help here.

 Thanks again for your help!
 Best,
 Debasish

[Aleksandra_K]

Could you also share how you compile/run the code? I tried your example and for me the algorithm converged as expected.

[DebasishBanerjee]

Oh, it is great to know that the code is running for you! Maybe I am making a mistake in the compilation. Here is my makefile

CXX=icpx
LIBS=-qmkl -lmkl_rt
OBJS=feast_ex2.o
OBJS1=feast_ex3.o

all: feast_ex2 feast_ex3

feast_ex2: ${OBJS}
${CXX} -o feast_ex2 ${OBJS} ${LIBS}

feast_ex3: ${OBJS1}
${CXX} -o feast_ex3 ${OBJS1} ${LIBS}
clean:
rm -rf feast_ex2 feast_ex3 ${OBJS} ${OBJS1} *.optrpt