bool My_function(std::vector &ev, Matrix &A, Matrix &B) { ev.resize(A->GetRowCount());//Set EigenValues massive size equal to matrix size //In my test it's equal to 8 MKL_INT ijob = -1; //Job indicator variable. On entry, a call to ?feast_?rci with ijob=-1 initializes the eigensolver. MKL_INT n = A->GetRowCount();//Sets the size of the problem. MKL_INT m0 = ev.size(); //On entry, specifies the initial guess for subspace dimension to be used. //m0 ? m where m is the total number of eigenvalues located in the interval [emin, emax]. //If the initial guess is wrong, Extended Eigensolver routines return info=3. MKL_INT m = ev.size(); //The total number of eigenvalues found in the interval [emin, emax]: 0 < m < m0. double epsout = 0.0; //On output, contains the relative error on the trace: |tracei - tracei-1| /max(|emin|, |emax|) double emin = 1e+4; //The lower bounds of the interval to be searched for eigenvalues; emin < emax. double emax = 2e+6; //The upper bounds of the interval to be searched for eigenvalues; emin < emax. MKL_INT info = 0; //On output, if info=0, the execution is successful. If info ? 0, see Output Eigensolver info Details. MKL_INT loop = 0; //On output, contains the number of refinement loop executed. Ignored on input. double *work; //Array of dimension n by m0, used to supply the temporary space for the ?feast_?rci computations. //Specifically, it contains the results of matrix-matrix multiplications. work = new double[n*m0]; double *aq; //Array of dimension n by m0, used to supply the temporary space for the ?feast_?rci computations. //Specifically, they contain the matrices for the reduced eigenvalue problem. aq = new double[n*m0]; double *sq; // --- || --- || --- sq = new double[n*m0]; double *q; //On entry, if fpm(5)=1, the array q(n, m) contains a basis of guess subspace where n is the order of the input matrix. q = new double[n*m]; double *res; //Array of length m0. On exit, the first m components contain the relative residual vector res = new double[m0]; MKL_INT fpm[128]; //see http://software.intel.com/sites/products/documentation/doclib/mkl_sa/11/mklman/index.htm feastinit (fpm); fpm[0] = 1; //Extended Eigensolver routines print runtime status to the screen. //fpm[1] = 8; //The number of contour points Ne = 8 (see the description of Feast algorithm). Must be one of {3,4,5,6,8,10,12,16,20,24,32,40,48}. //fpm[2] = 12; //Error trace double precision stopping criteria ? (? = 10-fpm(3)). //fpm[3] = 20; //Maximum number of Feast refinement loop allowed. If no convergence is reached within fpm(4) refinement loops, Extended Eigensolver routines return info=2. //fpm[4] = 0; //User initial subspace. If fpm(5)=0 then Extended Eigensolver routines generate initial subspace, if fpm(5)=1 the user supplied initial subspace is used. //fpm[5] = 0; //Feast stopping test. (see link above) //fpm[6] = 5; //Error trace single precision stopping criteria (10-fpm(7)) . //fpm[10] = 0; //Extended Eigensolver routines perform full contour integration. //fpm[13] = 0; //Return only subspace x after the first contour integration (0 or 1). //fpm[63] = 0; //Extended Eigensolver routines use the Intel MKL PARDISO default iparm settings defined by calling the pardisoinit subroutine. MKL_Complex16 ze; //Defines the coordinate along the complex contour. All values of ze are generated by ?feast_?rci internally. MKL_Complex16 *workc; //Array of dimension n by m0, used to supply the temporary space for the ?feast_?rci computations. //Specifically, it contains the solutions of various linear systems with complex coefficients. workc = new MKL_Complex16[m0]; cout << "\n---> phase " << ijob << " <---\n"; dfeast_srci(&ijob, &n, &ze, work, workc, aq, sq, fpm, &epsout, &loop, &emin, &emax, &m0, &ev[0], q, &m, res, &info); while (ijob !=0) { switch (ijob) { case 10 : //Factorize the complex matrix (ZeA-B) cout << "\n---> phase 10 <---\n"; break; case 11 : //Solve the complex linear system (ZeB-A)x=workc(1:N,1:M0) result in workc cout << "\n---> phase 11 <---\n"; break; case 30 : //Perform multiplication A*x(1:N,i:j) result in work(1:N,i:j) where i=fpm(25) and j=fpm(25)+fpm(24)?1 cout << "\n---> phase 30 <---\n"; break; case 40 : //Perform multiplication B*x(1:N,i:j) result in work(1:N,i:j) where i=fpm(25) and j=fpm(25)+fpm(24)?1 cout << "\n---> phase 40 <---\n"; break; default : //20 and 21 for complex problem only cout << "\n---> Something wrong! <---\n"; break; } dfeast_srci(&ijob, &n, &ze, work, workc, aq, sq, fpm, &epsout, &loop, &emin, &emax, &m0, &ev[0], q, &m, res, &info); } delete work; delete aq; delete sq; delete q; delete res; return true; }