#include "PardisoProvider.h"
#include "PerfCounter.h"
#include "cholmod.h"
#include <algorithm>
#include <iostream>
#include <mkl.h>

using namespace std;
using namespace Monitel::MathNative::Tools;
using namespace Monitel::MathNative::Pardiso;

Monitel::MathNative::Pardiso::PardisoProvider::PardisoProvider()
{	
	mkl_domain_set_num_threads(1, MKL_DOMAIN_PARDISO);
	ProviderInit();
}
/// <summary>
/// Инициализация решателя
/// </summary>
/// <param name="n">Число уравнений</param>
/// <param name="nrhs">Число правых частей</param>
/// <param name="a">Матрица</param>
/// <param name="ja">Индексы столбцов</param>
/// <param name="ia">Указатели строк</param>
/// <param name="executionTime">Время выполнения</param>
/// <returns></returns>
int Monitel::MathNative::Pardiso::PardisoProvider::Initialization(int n, int nrhs, Monitel::MathNative::Pardiso::Block2x2 * a, int * ja, int * ia, float *executionTime)
{	
	PerfCounter timer;
	double _ddum = 0;
	_phase = 11; //фаза инициализации
	_nrhs = nrhs;
	_mtype = MatrixType::RealNonsymmetric;
	ClearInternalMemory();
	UsingPermutationVector(2); //получаем вектор перестановок в фазе 1 - значение = 2

	timer.Start();
	_error = ConvertMatrix(n, a, ja, ia);
	if (_error == 0)
	{
		_permutations = new int[_internalRealMatrix.m];
		pardiso(_pt, &_maxfct, &_mnum, &_mtype, &_phase,
				&_internalRealMatrix.m, _internalRealMatrix.a, _internalRealMatrix.ia, _internalRealMatrix.ja,
				_permutations, &_nrhs, _iparm, &_msglvl, &_ddum, &_ddum, &_error);
	} 
	*executionTime = timer.Stop() * 1000;

	return _error;
}
int Monitel::MathNative::Pardiso::PardisoProvider::Initialization(int n, int nrhs, MKL_Complex16 * a, int * ja, int * ia, float * executionTime)
{
	PerfCounter timer;
	double _ddum = 0;
	int nnz = 0;
	_phase = 11; //фаза инициализации
	_nrhs = nrhs;
	_mtype = MatrixType::ComplexNonsymmetric;

	ClearInternalMemory();
	UsingPermutationVector(2); //получаем вектор перестановок в фазе 1 - значение = 2

	timer.Start();

	nnz = ia[n];
	_permutations = new int[n];

	_internalComplexMatrix.m = n;
	_internalComplexMatrix.nnz = nnz;
	_internalComplexMatrix.ac = new MKL_Complex16[nnz];
	_internalComplexMatrix.ja = new int[nnz];
	_internalComplexMatrix.ia = new int[n + 1];
	copy(a, a + nnz, _internalComplexMatrix.ac);
	copy(ja, ja + nnz, _internalComplexMatrix.ja);
	copy(ia, ia + (n + 1), _internalComplexMatrix.ia);

	pardiso(_pt, &_maxfct, &_mnum, &_mtype, &_phase,
			&_internalComplexMatrix.m, _internalComplexMatrix.ac, _internalComplexMatrix.ia, _internalComplexMatrix.ja,
			_permutations, &_nrhs, _iparm, &_msglvl, &_ddum, &_ddum, &_error);

	*executionTime = timer.Stop() * 1000;

	return _error;
}
int Monitel::MathNative::Pardiso::PardisoProvider::Initialization(int n, int nrhs, double * a, int * ja, int * ia, float * executionTime)
{
	PerfCounter timer;
	double _ddum = 0;
	int nnz = 0;

	_phase = 11; //фаза инициализации
	_nrhs = nrhs;
	_mtype = MatrixType::RealNonsymmetric;

	ClearInternalMemory();
	UsingPermutationVector(2); //получаем вектор перестановок в фазе 1 - значение = 2

	timer.Start();

	nnz = ia[n];
	_permutations = new int[n];

	_internalRealMatrix.m = n;
	_internalRealMatrix.nnz = nnz;
	_internalRealMatrix.a = new double[nnz];
	_internalRealMatrix.ja = new int[nnz];
	_internalRealMatrix.ia = new int[n + 1];
	copy(a, a + nnz, _internalRealMatrix.a);
	copy(ja, ja + nnz, _internalRealMatrix.ja);
	copy(ia, ia + (n + 1), _internalRealMatrix.ia);

	pardiso(_pt, &_maxfct, &_mnum, &_mtype, &_phase,
			&_internalRealMatrix.m, _internalRealMatrix.a, _internalRealMatrix.ia, _internalRealMatrix.ja,
			_permutations, &_nrhs, _iparm, &_msglvl, &_ddum, &_ddum, &_error);

	*executionTime = timer.Stop() * 1000;

	return _error;
}
/// <summary>
/// Факторизация матрицы
/// </summary>
/// <param name="PivotingPerturbation">Минимальное значение диагонального элемента</param>
/// <returns></returns>
int Monitel::MathNative::Pardiso::PardisoProvider::Factorization(int n, Monitel::MathNative::Pardiso::Block2x2 * a, int * ja, int * ia, 
																 float *converterTime, float *factorizationTime, int PivotingPerturbation)
{
	PerfCounter timer;
	double _ddum = 0;
	_phase = 22; //фаза факторизации
	ClearInternalRealMatrix();
	UsingPermutationVector(1); //используем полученный вектор перестановок в фазе 1 - значение = 1
	SetMinPivotValue(PivotingPerturbation);

	timer.Start();
	_error = ConvertMatrix(n, a, ja, ia);
	*converterTime = timer.Stop() * 1000;

	if (_error == 0)
	{
		timer.Start();
		pardiso(_pt, &_maxfct, &_mnum, &_mtype, &_phase,
				&_internalRealMatrix.m, _internalRealMatrix.a, _internalRealMatrix.ia, _internalRealMatrix.ja,
				_permutations, &_nrhs, _iparm, &_msglvl, &_ddum, &_ddum, &_error);
		*factorizationTime = timer.Stop() * 1000;
	}
	return _error;
}
int Monitel::MathNative::Pardiso::PardisoProvider::Factorization(int n, MKL_Complex16 * a, int * ja, int * ia, float * executionTime, int PivotingPerturbation)
{
	PerfCounter timer;
	double _ddum = 0;
	int nnz = 0;

	_phase = 22; //фаза факторизации
	ClearInternalComplexMatrix();
	UsingPermutationVector(1); //используем полученный вектор перестановок в фазе 1 - значение = 1
	SetMinPivotValue(PivotingPerturbation);

	nnz = ia[n];
	_internalComplexMatrix.m = n;
	_internalComplexMatrix.nnz = nnz;
	_internalComplexMatrix.ac = new MKL_Complex16[nnz];
	_internalComplexMatrix.ja = new int[nnz];
	_internalComplexMatrix.ia = new int[n + 1];
	copy(a, a + nnz, _internalComplexMatrix.ac);
	copy(ja, ja + nnz, _internalComplexMatrix.ja);
	copy(ia, ia + (n + 1), _internalComplexMatrix.ia);

	timer.Start();
	pardiso(_pt, &_maxfct, &_mnum, &_mtype, &_phase,
			&_internalComplexMatrix.m, _internalComplexMatrix.ac, _internalComplexMatrix.ia, _internalComplexMatrix.ja,
			_permutations, &_nrhs, _iparm, &_msglvl, &_ddum, &_ddum, &_error);
	*executionTime = timer.Stop() * 1000;

	return _error;
}
int Monitel::MathNative::Pardiso::PardisoProvider::Factorization(int n, double * a, int * ja, int * ia, float * executionTime, int PivotingPerturbation)
{
	PerfCounter timer;
	double _ddum = 0;
	int nnz = 0;

	_phase = 22; //фаза факторизации
	ClearInternalRealMatrix();
	UsingPermutationVector(1); //используем полученный вектор перестановок в фазе 1 - значение = 1
	SetMinPivotValue(PivotingPerturbation);

	nnz = ia[n];
	_internalRealMatrix.m = n;
	_internalRealMatrix.nnz = nnz;
	_internalRealMatrix.a = new double[nnz];
	_internalRealMatrix.ja = new int[nnz];
	_internalRealMatrix.ia = new int[n + 1];
	copy(a, a + nnz, _internalRealMatrix.a);
	copy(ja, ja + nnz, _internalRealMatrix.ja);
	copy(ia, ia + (n + 1), _internalRealMatrix.ia);

	timer.Start();
	pardiso(_pt, &_maxfct, &_mnum, &_mtype, &_phase,
			&_internalRealMatrix.m, _internalRealMatrix.a, _internalRealMatrix.ia, _internalRealMatrix.ja,
			_permutations, &_nrhs, _iparm, &_msglvl, &_ddum, &_ddum, &_error);
	*executionTime = timer.Stop() * 1000;

	return _error;
}
int Monitel::MathNative::Pardiso::PardisoProvider::Solve(MKL_Complex16 * rhs, MKL_Complex16 * solution, float * executionTime)
{
	PerfCounter timer;
	_phase = 33; //фаза решения СЛАУ
	UsingPermutationVector(1); //используем полученный вектор перестановок в фазе 1 - значение = 1

	timer.Start();
	pardiso(_pt, &_maxfct, &_mnum, &_mtype, &_phase,
			&_internalComplexMatrix.m, _internalComplexMatrix.ac, _internalComplexMatrix.ia, _internalComplexMatrix.ja,
			_permutations, &_nrhs, _iparm, &_msglvl, rhs, solution, &_error);
	*executionTime = timer.Stop() * 1000;

	return _error;
}
/// <summary>
/// Решение СЛАУ
/// </summary>
/// <param name="rhs">Правая часть системы уравнений</param>
/// <param name="solution">Результат решения СЛАУ</param>
/// <param name="executionTime"></param>
/// <returns></returns>
int Monitel::MathNative::Pardiso::PardisoProvider::Solve(double * rhs, double * solution, float *executionTime)
{
	PerfCounter timer;
	_phase = 33; //фаза решения СЛАУ
	UsingPermutationVector(1); //используем полученный вектор перестановок в фазе 1 - значение = 1

	timer.Start();
	pardiso(_pt, &_maxfct, &_mnum, &_mtype, &_phase,
			&_internalRealMatrix.m, _internalRealMatrix.a, _internalRealMatrix.ia, _internalRealMatrix.ja,
			_permutations, &_nrhs, _iparm, &_msglvl, rhs, solution, &_error);
	*executionTime = timer.Stop() * 1000;
	
	return _error;
}
/// <summary>
/// Установить метод перестановок
/// </summary>
/// <param name="method"></param>
void Monitel::MathNative::Pardiso::PardisoProvider::SetPermutationMethod(int method)
{
	if (method >= 0 & method <= 3)
	{
		_iparm1 = method;
	}	
	else
		_iparm1 = 2;

	SetInitParam();
}
/// <summary>
/// Установить число шагов итеративного решателя
/// </summary>
/// <param name="steps">Число шагов</param>
void Monitel::MathNative::Pardiso::PardisoProvider::SetIterativeRefinementStepNumber(int steps)
{
	_iparm7 = steps;
	SetInitParam();
}
void Monitel::MathNative::Pardiso::PardisoProvider::SetUseAVX(bool status)
{
	_useAVX = status;

	if (!_useAVX)
	{
		mkl_cbwr_set(MKL_CBWR_SSSE3);
	}
}
Monitel::MathNative::Pardiso::PardisoProvider::~PardisoProvider()
{
	double ddum = 0;
	_phase = -1;           /* Release internal memory. */
	if (_mtype == MatrixType::RealNonsymmetric)
	{
		pardiso(_pt, &_maxfct, &_mnum, &_mtype, &_phase,
				&_internalRealMatrix.m, &ddum, _internalRealMatrix.ia, _internalRealMatrix.ja, _permutations, &_nrhs,
				_iparm, &_msglvl, &ddum, &ddum, &_error);
	} else if (_mtype == MatrixType::ComplexNonsymmetric)
	{
		pardiso(_pt, &_maxfct, &_mnum, &_mtype, &_phase,
				&_internalComplexMatrix.m, &ddum, _internalComplexMatrix.ia, _internalComplexMatrix.ja, _permutations, &_nrhs,
				_iparm, &_msglvl, &ddum, &ddum, &_error);
	}

	ClearInternalMemory();
}
/// <summary>
/// Инициализация провайдера
/// </summary>
void Monitel::MathNative::Pardiso::PardisoProvider::ProviderInit()
{
	for (size_t i = 0; i < 64; i++)
	{
		_iparm[i] = 0;
	}

	SetInitParam();

	for (size_t i = 0; i < 64; i++)
	{
		_pt[i] = 0;
	}
}
/// <summary>
/// Установить значения в вектор инициализации
/// </summary>
void Monitel::MathNative::Pardiso::PardisoProvider::SetInitParam()
{
	_iparm[0] = 1;         /* No solver default */
	_iparm[1] = _iparm1;
	_iparm[3] = _iparm3;
	_iparm[4] = _iparm4;
	_iparm[7] = _iparm7;
	_iparm[9] = _iparm9;
	_iparm[10] = _iparm10;
	_iparm[12] = _iparm12;
	_iparm[20] = _iparm20;
	_iparm[23] = _iparm23;
	_iparm[33] = _iparm33;
	_iparm[34] = _iparm34;
}
/// <summary>
/// Использование вектора перестановок
/// </summary>
/// <param name="mode"></param>
void Monitel::MathNative::Pardiso::PardisoProvider::UsingPermutationVector(int mode)
{
	if (mode >= 0 & mode <= 2)
	{
		_iparm4 = mode;
	} else
		_iparm4 = 2;

	SetInitParam();
}
void Monitel::MathNative::Pardiso::PardisoProvider::SetMinPivotValue(int level)
{
	_iparm9 = level;
	SetInitParam();
}
/// <summary>
/// Очистить внутренюю память
/// </summary>
void Monitel::MathNative::Pardiso::PardisoProvider::ClearInternalMemory()
{
	ClearInternalRealMatrix();
	ClearInternalComplexMatrix();

	if (_permutations != nullptr)
		delete[] _permutations;
	_permutations = nullptr;
}
void Monitel::MathNative::Pardiso::PardisoProvider::ClearInternalRealMatrix()
{
	if (_internalRealMatrix.a != nullptr)
		delete[] _internalRealMatrix.a;
	if (_internalRealMatrix.ja != nullptr)
		delete[] _internalRealMatrix.ja;
	if (_internalRealMatrix.ia != nullptr)
		delete[] _internalRealMatrix.ia;

	_internalRealMatrix.a = nullptr;
	_internalRealMatrix.ja = nullptr;
	_internalRealMatrix.ia = nullptr;
}
void Monitel::MathNative::Pardiso::PardisoProvider::ClearInternalComplexMatrix()
{
	if (_internalComplexMatrix.ac != nullptr)
		delete[] _internalComplexMatrix.a;
	if (_internalComplexMatrix.ja != nullptr)
		delete[] _internalComplexMatrix.ja;
	if (_internalComplexMatrix.ia != nullptr)
		delete[] _internalComplexMatrix.ia;

	_internalComplexMatrix.ac = nullptr;
	_internalComplexMatrix.ja = nullptr;
	_internalComplexMatrix.ia = nullptr;
}
/// <summary>
/// Конвертировать блочную матрицу в плоскую
/// </summary>
/// <param name="n">Число уравнений</param>
/// <param name="matrix">Матрица</param>
/// <param name="ja">Индексы столбцов</param>
/// <param name="ia">Указатели строк</param>
/// <param name="outn">Выходное число уравнений</param>
/// <param name="outMatrix">Выходная матрица</param>
/// <param name="outja">Выходные индексы столбцов</param>
/// <param name="outia">Выходные указатели строк</param>
int Monitel::MathNative::Pardiso::PardisoProvider::ConvertMatrix(MKL_INT n, Block2x2 * matrix, MKL_INT * ja, MKL_INT * ia)
{
	try
	{
		cholmod_sparse *A = 0;
		cholmod_sparse *B = 0;
		cholmod_sparse *AB_tmp;	//промежуточная матрица
		cholmod_sparse *C = 0;
		cholmod_sparse *D = 0;
		cholmod_sparse *CD_tmp; //промежуточная матрица
		cholmod_sparse *J;		//Итоговая матрица

		cholmod_common c;
		int values = 1;
		int _M = n;				//размер матрицы М*М
		int nnzmax = ia[n];	//число ненулевых элементов в матрице

		double *csrMatrix_A = nullptr;
		double *csrMatrix_B = nullptr;
		double *csrMatrix_C = nullptr;
		double *csrMatrix_D = nullptr;

		int *csrRowPtr;
		int *csrColInd;

		double *cscMatrix_A = nullptr;
		double *cscMatrix_B = nullptr;
		double *cscMatrix_C = nullptr;
		double *cscMatrix_D = nullptr;

		int *cscRowInd = nullptr;
		int *cscColPtr = nullptr;
		int *job = nullptr;

		A = new cholmod_sparse;
		B = new cholmod_sparse;
		C = new cholmod_sparse;
		D = new cholmod_sparse;

		csrMatrix_A = new double[nnzmax];
		csrMatrix_B = new double[nnzmax];
		csrMatrix_C = new double[nnzmax];
		csrMatrix_D = new double[nnzmax];

		csrRowPtr = ia;
		csrColInd = ja;

		cscMatrix_A = new double[nnzmax];
		cscMatrix_B = new double[nnzmax];
		cscMatrix_C = new double[nnzmax];
		cscMatrix_D = new double[nnzmax];

		cscRowInd = new int[nnzmax];
		cscColPtr = new int[_M + 1];

		//растасовываем структуру по матрицам
		for (int i = 0; i < nnzmax; i++)
		{
			csrMatrix_A[i] = matrix[i].f00;
			csrMatrix_B[i] = matrix[i].f01;
			csrMatrix_C[i] = matrix[i].f10;
			csrMatrix_D[i] = matrix[i].f11;
		}

		//Транспонируем матрицы используя MKL. Cholmod работает с форматом CSC.	
		job = new int[8];
		int info;
		fill(job, job + 8, 0);

		job[1 - 1] = 0; //CSR to CSC
		job[2 - 1] = 0; //нумерация с нуля для CSR
		job[3 - 1] = 0; //нумерация с нуля для СSC
		job[6 - 1] = 1; //заполняем все массивы

		mkl_dcsrcsc(job, &_M, csrMatrix_A, csrColInd, csrRowPtr, cscMatrix_A, cscRowInd, cscColPtr, &info);
		mkl_dcsrcsc(job, &_M, csrMatrix_B, csrColInd, csrRowPtr, cscMatrix_B, cscRowInd, cscColPtr, &info);
		mkl_dcsrcsc(job, &_M, csrMatrix_C, csrColInd, csrRowPtr, cscMatrix_C, cscRowInd, cscColPtr, &info);
		mkl_dcsrcsc(job, &_M, csrMatrix_D, csrColInd, csrRowPtr, cscMatrix_D, cscRowInd, cscColPtr, &info);

		cholmod_start(&c);			   //  start CHOLMOD

		//Заполняем А
		A->nrow = _M;
		A->ncol = _M;
		A->nzmax = nnzmax;
		A->p = cscColPtr;
		A->i = cscRowInd;
		A->x = cscMatrix_A;
		A->stype = 0;
		A->itype = CHOLMOD_INTLONG;
		A->xtype = CHOLMOD_REAL;
		A->dtype = CHOLMOD_DOUBLE;
		A->sorted = true;
		A->packed = true;

		//Заполняем В
		B->nrow = _M;
		B->ncol = _M;
		B->nzmax = nnzmax;
		B->p = cscColPtr;
		B->i = cscRowInd;
		B->x = cscMatrix_B;
		B->stype = 0;
		B->itype = CHOLMOD_INTLONG;
		B->xtype = CHOLMOD_REAL;
		B->dtype = CHOLMOD_DOUBLE;
		B->sorted = true;
		B->packed = true;

		//Заполняем С
		C->nrow = _M;
		C->ncol = _M;
		C->nzmax = nnzmax;
		C->p = cscColPtr;
		C->i = cscRowInd;
		C->x = cscMatrix_C;
		C->stype = 0;
		C->itype = CHOLMOD_INTLONG;
		C->xtype = CHOLMOD_REAL;
		C->dtype = CHOLMOD_DOUBLE;
		C->sorted = true;
		C->packed = true;

		//Заполняем D
		D->nrow = _M;
		D->ncol = _M;
		D->nzmax = nnzmax;
		D->p = cscColPtr;
		D->i = cscRowInd;
		D->x = cscMatrix_D;
		D->stype = 0;
		D->itype = CHOLMOD_INTLONG;
		D->xtype = CHOLMOD_REAL;
		D->dtype = CHOLMOD_DOUBLE;
		D->sorted = true;
		D->packed = true;

		AB_tmp = cholmod_horzcat(A, B, values, &c);  // AB_tmp = A|B
		CD_tmp = cholmod_horzcat(C, D, values, &c);  // CD_tmp = C|D
		J = cholmod_vertcat(AB_tmp, CD_tmp, values, &c); //итоговая матрица

		//размерность итоговой матрицы
		_internalRealMatrix.m = J->nrow;
		//число ненулевых элементов в итоговой матрице
		_internalRealMatrix.nnz = J->nzmax;
		//Память по итоговую матрицу
		_internalRealMatrix.a = new double[J->nzmax];
		_internalRealMatrix.ja = new MKL_INT[J->nzmax];
		_internalRealMatrix.ia = new MKL_INT[J->nrow + 1];

		fill(job, job + 8, 0);

		job[1 - 1] = 1; //CSC to CSR
		job[2 - 1] = 0; //нумерация с нуля для CSR
		job[3 - 1] = 0; //нумерация с нуля для СSC
		job[6 - 1] = 1; //заполняем все массивы

		mkl_dcsrcsc(job,
					&_internalRealMatrix.m, _internalRealMatrix.a, _internalRealMatrix.ja, _internalRealMatrix.ia,
					(double*)J->x, (int*)J->i, (int*)J->p, &info);

		cholmod_free_sparse(&AB_tmp, &c);
		cholmod_free_sparse(&CD_tmp, &c);
		cholmod_free_sparse(&J, &c);

		free(csrMatrix_A);
		free(csrMatrix_B);
		free(csrMatrix_C);
		free(csrMatrix_D);

		free(A);
		free(B);
		free(C);
		free(D);

		free(cscMatrix_A);
		free(cscMatrix_B);
		free(cscMatrix_C);
		free(cscMatrix_D);

		free(cscRowInd);
		free(cscColPtr);

		free(job);
		cholmod_finish(&c);			    // finish CHOLMOD

		return 0;
	}
	catch (const std::exception& ex)
	{
		return -1;
	}
}
int Monitel::MathNative::Pardiso::PardisoProvider::GetCBWRBranch()
{
	int my_cbwr_branch;

	my_cbwr_branch = mkl_cbwr_get(MKL_CBWR_BRANCH);

	return my_cbwr_branch;
}
int * Monitel::MathNative::Pardiso::PardisoProvider::GetPermutations()
{
	return _permutations;
}
void * Monitel::MathNative::Pardiso::PardisoProvider::GetFactorizedMatrixDiagonal()
{
	MKL_INT mnum = 1;
	MKL_INT error = 0;

	void* df = nullptr;
	void* da = nullptr;

	pardiso_getdiag(_pt, df,da, &mnum, &error);
	if (error != 0)
	{
		return df;
	} else
	{
		return nullptr;
	}	
}
void * Monitel::MathNative::Pardiso::PardisoProvider::GetInitialMatrixDiagonal()
{
	MKL_INT mnum = 1;
	MKL_INT error = 0;

	void* df = nullptr;
	void* da = nullptr;

	pardiso_getdiag(_pt, df, da, &mnum, &error);
	if (error != 0)
	{
		return da;
	} else
	{
		return nullptr;
	}
}
INativePardisoProvider* Monitel::MathNative::Pardiso::CreatePardisoSolverInstanse()
{
	INativePardisoProvider *newIns = new PardisoProvider();
	return newIns;
	//return nullptr;
}
void Monitel::MathNative::Pardiso::ReleasePardisoSolverInstanse(INativePardisoProvider *pardisoProvider)
{
	delete pardisoProvider;
	pardisoProvider = nullptr;
}

