bool
My_function(std::vector<double> &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;
}