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