#pragma once

#include <iostream>


#include <mkl_rci.h>
#include <mkl_blas.h>
#include <mkl_spblas.h>
#include <mkl_service.h>
#include <stdio.h>


using namespace std;


template <class T>
class SPMV  {

protected:// shared variable amoung children
	
	void call_gmres_solver_mkl(const vector<T>& coefNeigbor, const vector<T>& RHS,
		vector<T>& phi, vector<T>& phiOld);

};




template<class T>
void SPMV<T>::call_gmres_solver_mkl(const vector<T>& coefNeigbor, const vector<T>& RHS, 
	vector<T>& phi, vector<T>& phiOld) {

	if (sizeof(T) == sizeof(float)) {
		cout << "currently this program is written for double floatting point numbers only.\nplease update the source code and recompiler" << endl;
		cout << "Program will terminate now" << endl;
		cin.get();
		exit(0);
	}
	

	MKL_INT n = nFl;// number of unknows
	MKL_INT nEle = (nFl + (nFl*nface) - (nTtl - nFl));// digonal elements + #of internal faces

	if (nEle <= 0) {
		cout << "The number of interfaces must be greater than zero. current value =" << n << endl;
		cout << "Program will exit now" << endl;
		cin.get();
		exit(0);
	}

	MKL_INT RCI_request, itercount, expected_itercount = 8;// , i;

	/* Fill all arrays containing matrix data. */
	MKL_INT *ia = new MKL_INT[n + 1];
	MKL_INT *ja = new MKL_INT[nEle];


	// create matrix A
	int id = 0;
	ia[0] = 0;

	int n0 = n * (2 * n + 1) + (n * (n + 9)) / 2 + 1;
	double *tmp = (double *)malloc(n0 * sizeof(double));
	double *rhs = (double *)malloc(n * sizeof(double));
	double *solution = (double *)malloc(n * sizeof(double));

	MKL_INT ipar[128];
	double dpar[128];

	fill_n(tmp, n0, 0.0);
	
	for (int i = 0; i < n; i++) {
		solution[i] = phiOld[i];// initial guess is the old values
	}


	/* Initialize the initial guess*/
	for (int i = 0; i < n; i++) {
		solution[i] = phiOld[i];
		rhs[i] = RHS[i];
	}
	
	/*------------------- Initialize the solver---------------------------------*/
	dfgmres_init(&n, solution, rhs, &RCI_request, ipar, dpar, tmp);
	if (RCI_request != 0) {
		cout << "Some error occured during the gmres initialisation phase. Error code =" << RCI_request << endl;
		goto FAILED;
	}
	
	/*-------------------------------------------------------------------------*/
FAILED:printf("\nThis example FAILED as the solver has returned the ERROR code %d", RCI_request);



}