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;
}