#pragma once

#include <mkl.h>
#include <mkl_types.h>
#include "IPardisoProvider.h"

namespace Monitel
{
	namespace MathNative
	{
		namespace Pardiso
		{
			/// <summary>
			/// Провайдер Pardiso
			/// </summary>
			class PARDISO_API PardisoProvider : public INativePardisoProvider
			{
				/// <summary>
				/// csr матрица
				/// </summary>
				struct csrMatrix
				{
					/// <summary>
					/// Число строк матрицы
					/// </summary>
					MKL_INT m = 0;
					/// <summary>
					/// Число ненулевых элементов
					/// </summary>
					MKL_INT nnz = 0;
					/// <summary>
					/// Матрица
					/// </summary>
					double *a = nullptr;
					/// <summary>
					/// Матрица
					/// </summary>
					MKL_Complex16 *ac = nullptr;
					/// <summary>
					/// Индексы столбцов
					/// </summary>
					MKL_INT *ja = nullptr;
					/// <summary>
					/// Указатели строк
					/// </summary>
					MKL_INT *ia = nullptr;
				};
				/// <summary>
				/// Тип матрицы
				/// </summary>
				enum MatrixType
				{
					RealNonsymmetric = 11,
					ComplexNonsymmetric = 13
				};

			public:
				/// <summary>
				/// Конструктор
				/// </summary>
				PardisoProvider();
				/// <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 Initialization(int n, int nrhs, Monitel::MathNative::Pardiso::Block2x2 * a, int * ja, int * ia, float *executionTime);
				/// <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 Initialization(int n, int nrhs, MKL_Complex16 * a, int * ja, int * ia, float *executionTime);
				/// <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 Initialization(int n, int nrhs, double * a, int * ja, int * ia, float *executionTime);
				/// <summary>
				/// Факторизация
				/// </summary>
				/// <param name="n">Число уравнений</param>
				/// <param name="a">Матрица</param>
				/// <param name="ja">Индексы столбцов</param>
				/// <param name="ia">Указатели строк</param>
				/// <param name="converterTime">Время конвертации матрицы</param>
				/// <param name="factorizationTime">Время факторизации матрицы</param>
				/// <param name="PivotingPerturbation">Минимальное значение диагонального элемента</param>
				/// <returns></returns>
				int Factorization(int n, Monitel::MathNative::Pardiso::Block2x2 * a, int * ja, int * ia, float *converterTime, float *factorizationTime, int PivotingPerturbation = 13);
				/// <summary>
				/// Факторизация
				/// </summary>
				/// <param name="n">Число уравнений</param>
				/// <param name="a">Матрица</param>
				/// <param name="ja">Индексы столбцов</param>
				/// <param name="ia">Указатели строк</param>
				/// <param name="executionTime">Время выполнения</param>
				/// <param name="PivotingPerturbation">Минимальное значение диагонального элемента</param>
				/// <returns></returns>
				int Factorization(int n, MKL_Complex16 * a, int * ja, int * ia, float *executionTime, int PivotingPerturbation = 13);
				/// <summary>
				/// Факторизация
				/// </summary>
				/// <param name="n">Число уравнений</param>
				/// <param name="a">Матрица</param>
				/// <param name="ja">Индексы столбцов</param>
				/// <param name="ia">Указатели строк</param>
				/// <param name="executionTime">Время выполнения</param>
				/// <param name="PivotingPerturbation">Минимальное значение диагонального элемента</param>
				/// <returns></returns>
				int Factorization(int n, double * a, int * ja, int * ia, float *executionTime, int PivotingPerturbation = 13);
				/// <summary>
				/// Решение СЛАУ
				/// </summary>
				/// <param name="rhs">Правая часть системы уравнений</param>
				/// <param name="solution">Результат решения СЛАУ</param>
				/// <param name="executionTime">Время выполнения</param>
				/// <returns></returns>
				int Solve(MKL_Complex16 * rhs, MKL_Complex16 * solution, float *executionTime);
				/// <summary>
				/// Решение СЛАУ
				/// </summary>
				/// <param name="rhs">Правая часть системы уравнений</param>
				/// <param name="solution">Результат решения СЛАУ</param>
				/// <param name="executionTime">Время выполнения</param>
				/// <returns></returns>
				int Solve(double * rhs, double * solution, float *executionTime);
				/// <summary>
				/// Получить вектор перестановок. Длина вектора равна n.
				/// </summary>
				/// <returns></returns>
				int * GetPermutations();
				/// <summary>
				/// Получить диагональ факторизованной матрицы. Длина вектора равна n.
				/// </summary>
				/// <returns></returns>
				void * GetFactorizedMatrixDiagonal();
				/// <summary>
				/// Получить диагональ исходной матрицы. Длина вектора равна n.
				/// </summary>
				/// <returns></returns>
				void * GetInitialMatrixDiagonal();
				/// <summary>
				/// Установить метод перестановок
				/// </summary>
				/// <param name="method"></param>
				void SetPermutationMethod(int method);
				/// <summary>
				/// Установить число шагов итеративного решателя
				/// </summary>
				void SetIterativeRefinementStepNumber(int steps);
				/// <summary>
				/// Установить использование AVX инструкций
				/// </summary>
				/// <param name="status"></param>
				void SetUseAVX(bool status);
				/// <summary>
				/// Деструктор
				/// </summary>
				~PardisoProvider();

			private:
				/// <summary>
				/// Инициализация провайдера
				/// </summary>
				void ProviderInit();
				/// <summary>
				/// Установить значения в вектор инициализации
				/// </summary>
				void SetInitParam();
				/// <summary>
				/// Использование вектора перестановок
				/// </summary>
				/// <param name="mode"></param>
				void UsingPermutationVector(int mode);
				/// <summary>
				/// Установить минимальное значение диагонального элемента
				/// </summary>
				/// <param name="level"></param>
				void SetMinPivotValue(int level);
				/// <summary>
				/// Очистить внутренюю память
				/// </summary>
				void ClearInternalMemory();
				/// <summary>
				/// Очистить внутренюю вещественную матрицу
				/// </summary>
				void ClearInternalRealMatrix();
				/// <summary>
				/// Очистить внутренюю комплексную матрицу
				/// </summary>
				void ClearInternalComplexMatrix();
				/// <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 ConvertMatrix(
					MKL_INT n,
					Monitel::MathNative::Pardiso::Block2x2 *matrix, MKL_INT *ja, MKL_INT *ia);
				int GetCBWRBranch();

				/// <summary>
				/// Тип матрицы
				/// 11* - вещественная несемметричная
				/// 13 - комплексная несемметричная
				/// </summary>
				MKL_INT _mtype = 11;
				/// <summary>
				/// Число правых частей
				/// </summary>
				MKL_INT _nrhs = 1;
				/// <summary>
				/// Указатель на внутренюю память
				/// </summary>
				void *_pt[64];
				/// <summary>
				/// Настроечные параметры PARDISO
				/// </summary>
				MKL_INT _iparm[64];
				/// <summary>
				/// Число наборов факторов
				/// </summary>
				MKL_INT _maxfct = 1;
				/// <summary>
				/// Номер набора для факторизации
				/// </summary>
				MKL_INT _mnum = 1;
				/// <summary>
				/// Уровень отладки
				/// 0 - ничего не печатать в консоль
				/// 1 - печатать статистическую информацию в консоль
				/// </summary>
				MKL_INT _msglvl = 0;
				/// <summary>
				/// Флаг ошибки
				/// </summary>
				MKL_INT _error = 0;
				/// <summary>
				/// Фаза решателя
				/// 11* - инициализация
				/// 22 - факторизация
				/// 33 - решение СЛАУ
				/// </summary>
				MKL_INT _phase = 11;
				/// <summary>
				/// Число уравнений СЛАУ
				/// </summary>
				MKL_INT _n = 0;
				/// <summary>
				/// Выбор метода перестановок
				/// 0 - The minimum degree algorithm
				/// 2* - The nested dissection algorithm from the METIS package
				/// 3 - The parallel (OpenMP) version of the nested dissection algorithm. 
				/// It can decrease the time of computations on multi-core computers, 
				/// especially when Intel MKL PARDISO Phase 1 takes significant time.
				/// </summary>
				MKL_INT _iparm1 = 2;
				/// <summary>
				/// Preconditioned CGS/CG
				/// </summary>
				MKL_INT _iparm3 = 0;
				/// <summary>
				/// Использование перестановок
				/// 0 - User permutation in the perm array is ignored.
				/// 1 - Intel MKL PARDISO uses the user supplied fill-in reducing permutation from the perm array. iparm[1] is ignored.
				/// 2* - Intel MKL PARDISO returns the permutation vector computed at phase 1 in the perm array.
				/// </summary>
				MKL_INT _iparm4 = 2;
				/// <summary>
				/// Iterative refinement step
				/// </summary>
				MKL_INT _iparm7 = -1;
				/// <summary>
				/// Pivoting perturbation
				/// </summary>
				MKL_INT _iparm9 = 13;
				/// <summary>
				/// Scaling vectors
				/// 0 - Disable scaling. Default for symmetric indefinite matrices.
				/// 1* - Enable scaling. Default for nonsymmetric matrices.
				/// </summary>
				MKL_INT _iparm10 = 1;
				/// <summary>
				/// Improved accuracy using (non-) symmetric weighted matching.
				/// 0 - Disable matching. Default for symmetric indefinite matrices.
				/// 1* - Enable matching. Default for nonsymmetric matrices.
				/// </summary>
				MKL_INT _iparm12 = 1;
				/// <summary>
				/// Pivoting for symmetric indefinite matrices.
				/// </summary>
				MKL_INT _iparm20 = 2;
				/// <summary>
				/// Parallel factorization control.
				/// 0* - Intel MKL PARDISO uses the classic algorithm for factorization.
				/// 1 - Intel MKL PARDISO uses a two-level factorization algorithm. 
				/// This algorithm generally improves scalability in case of parallel factorization on many OpenMP threads (more than eight).
				/// 10 - Intel MKL PARDISO uses an improved two-level factorization algorithm for nonsymmetric matrices.
				/// </summary
				MKL_INT _iparm23 = 0;
				/// <summary>
				/// Optimal number of OpenMP threads for conditional numerical reproducibility (CNR) mode
				/// </summary>
				MKL_INT _iparm33 = 1;
				/// <summary>
				/// Нумерация индексов столбцов и строк
				/// 0 - начинаются с единици (стиль языка Фортран)
				/// 1* - начинаются с нуля (стиль языка С)
				/// </summary>
				MKL_INT _iparm34 = 1;
				/// <summary>
				/// Использовать инструкции AVX
				/// </summary>
				bool _useAVX = false;

				/***************************************/
				/// <summary>
				/// Внутреняя вещественная матрица матрица
				/// </summary>
				csrMatrix _internalRealMatrix;
				/// <summary>
				/// Внутреняя комплексная матрица
				/// </summary>
				csrMatrix _internalComplexMatrix;
				/// <summary>
				/// Вектор перестановок
				/// </summary>
				MKL_INT *_permutations = nullptr;
			};
		}
	}
}


