#include <random>
#include <chrono>
#include <iostream>
#include <vector>
#include <functional>
#include <numeric>
#include <intrin.h>
#include <omp.h>
#include <iomanip>
using namespace std;
void printMatrix(int n, double* m);

void luDecomposition(int n, double* mat)
{
	double **lower = 0;
	double **upper = 0;
	lower = upper = new double *[n];
	for (int i = 0; i < n; i++)
	{
		lower[i] = new double[n];
		upper[i] = new double[n];
	}
	//double[n][n], upper[n][n];
	//memset(lower, 0, sizeof(lower));
	//memset(upper, 0, sizeof(upper));
	for (int i = 0; i < n; i++) {

		for (int k = i; k < n; k++) {

			double sum = 0;
			for (int j = 0; j < i; j++)
				sum += (lower[i][j] * upper[j][k]);

			upper[i][k] = mat[(n*i) + k] - sum;
		}

		for (int k = i; k < n; k++) {
			if (i == k)
				lower[i][i] = 1;
			else {

				double sum = 0;
				for (int j = 0; j < i; j++)
					sum += (lower[k][j] * upper[j][i]);

				lower[k][i] = (mat[(n*k) + i] - sum) / upper[i][i];
			}
		}
	}

	cout << setw(6) << "Lower Triangular" << endl << endl;


	for (int i = 0; i < n; i++) {
		for (int j = 0; j < n; j++)
		{
			cout << lower[i][j] << setw(6) << " ";
		}
		cout << endl;
	}
	cout << endl << endl;
	cout << setw(6) << "upper Triangular" << endl << endl;
	for (int i = 0; i < n; i++) {
		for (int j = 0; j < n; j++)
		{
			cout << upper[i][j] << setw(6) << " ";
		}
		cout << endl;
	}
	cout << endl << endl;

}

void lup_od_omp(double* a, int n) {

	int i, j, k;

	for (k = 0; k < n - 1; ++k)
	{
#pragma omp parallel for shared(a,n,k) private(i) schedule(static, 64)
		for (i = k + 1; i < n; i++)
		{
			a[i*n + k] /= a[k*n + k];
			for (j = k + 1; j < n; j++)
			{
				a[i*n + j] -= a[i*n + k] * a[k*n + j];
			}
		}
	}
}


void printMatrix(int n, double* m)
{
	for (size_t i = 0; i < n; ++i)
	{
		for (size_t j = 0; j < n; ++j)
		{
			cout.precision(4);
			std::cout << setw(6) << m[(n * i) + j] << setw(8) << " ";
		}
		std::cout << std::endl;
	}
	cout << endl;
}

int main()
{
	size_t runCount = 1000;
	std::vector<std::vector<long long>> result;
	result.emplace_back();
	result.emplace_back();
	result.emplace_back();

	size_t n = 10;

	double *m1 = nullptr;
	double *temp1 = nullptr;
	double *temp2 = nullptr;


	try
	{
		m1 = new double[n * n];
		temp1 = new double[n * n];
		temp2 = new double[n * n];
	}
	catch (std::bad_alloc& e)
	{
		std::cout << "Need more memory" << std::endl;
		delete[] m1;
		delete[] temp1;
		delete[] temp2;
		return 1;
	}

	for (size_t currentRun = 0; currentRun < runCount; ++currentRun)
	{
		std::cout << "current run:  " << currentRun << std::endl;
		std::random_device rd;

		for (size_t i = 0; i < n * n; ++i)
		{
			m1[i] = static_cast<double>(rd()) / static_cast<double>(rd());
		}
		std::cout << endl << "Matrix:\r\n";
		printMatrix(n, m1);

		std::chrono::high_resolution_clock::time_point t1 = std::chrono::high_resolution_clock::now();
		luDecomposition(n, m1);
		std::chrono::high_resolution_clock::time_point t2 = std::chrono::high_resolution_clock::now();

		auto duration = std::chrono::duration_cast<std::chrono::microseconds>(t2 - t1).count();

		result[0].push_back(duration);
		std::cout << "Seq duration: " << duration << std::endl;

		for (size_t i = 0; i < n * n; ++i)
		{
			temp1[i] = 0;
			temp2[i] = 0;
		}

		t1 = std::chrono::high_resolution_clock::now();
		lup_od_omp(m1, n);
		t2 = std::chrono::high_resolution_clock::now();

		duration = std::chrono::duration_cast<std::chrono::microseconds>(t2 - t1).count();

		result[1].push_back(duration);

		//for (size_t i = 0; i < n * n; ++i)
		//{
			//[i] = 0;
		//}

		//t1 = std::chrono::high_resolution_clock::now();
		//sumOpenMP(n * n, m1, m2, res);
		//t2 = std::chrono::high_resolution_clock::now();

		duration = std::chrono::duration_cast<std::chrono::microseconds>(t2 - t1).count();

		result[2].push_back(duration);
		std::cout << "Vec duration: " << duration << std::endl;
	}

	for (size_t i = 0; i < result.size(); ++i)
	{
		auto r = std::accumulate(result[i].begin(), result[i].end(), 0);
		std::cout << static_cast<double>(r) / static_cast<double>(result[i].size()) << std::endl;
	}

	//sumOpenMP(n * n, m1, m2, res);

	//std::cout << "m1:\r\n";
	//printMatrix(n, m1);

	//std::cout << "m2:\r\n";
	//printMatrix(n, m2);


	delete[] m1;
	delete[] temp1;
	delete[] temp2;


	std::cin >> n;

	return 0;
}