#include "pch.h"
#include "PardisoProvider.h"
#include <algorithm>
#include <iostream>
#include <mkl.h>

using namespace std;

PardisoProvider::PardisoProvider()
{
	mkl_domain_set_num_threads(1, MKL_DOMAIN_PARDISO);
	ProviderInit();
}

int PardisoProvider::Initialization(int n, int nrhs, double * a, int * ja, int * ia)
{
	double _ddum = 0;
	_nnz = 0;

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

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

	_n = n;
	_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);

	return _error;
}

int PardisoProvider::Factorization()
{
	double _ddum = 0;
	int nnz = 0;

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

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

	return _error;
}

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

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

	return _error;
}

PardisoProvider::~PardisoProvider()
{
	double ddum = 0;
	_phase = -1;           /* Release internal memory. */

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


	ClearInternalMemory();
}
/// <summary>
/// Инициализация провайдера
/// </summary>
void 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 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 PardisoProvider::UsingPermutationVector(int mode)
{
	if (mode >= 0 & mode <= 2)
	{
		_iparm4 = mode;
	} else
		_iparm4 = 2;

	SetInitParam();
}
/// <summary>
/// Очистить внутренюю память
/// </summary>
void PardisoProvider::ClearInternalMemory()
{
	ClearInternalRealMatrix();

	if (_permutations != nullptr)
		delete[] _permutations;
	_permutations = nullptr;
}
void 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;
}