#include "app.h"

#define SEQ_OUTPUT_FILENAME "out_seq.csv"
#define SEQ_PER_BLOCK_OUTPUT_FILENAME "out_seq_per_block.csv"
#define PARALLEL_OUTPUT_FILENAME "out_parallel.csv"
#define PARALLEL_PER_BLOCK_OUTPUT_FILENAME "out_parallel_per_block.csv"
#define SCENARIO_OUTPUT_FILENAME "out.csv"
using namespace std;
using namespace std::chrono;

int main(int argc, char **argv)
{

	//CreateData(SCENARIO_OUTPUT_FILENAME);
	scenario();


	/*int size, block_size, num_iter,max_threads;
	high_resolution_clock::time_point t1, t2;
	
	//saisie des parametres
	cin >> size >> block_size >> num_iter>>max_threads;
	//initialisation des matrices
	float **A = new float*[size];
	float **B = new float*[size];
	float **C = new float*[size];
	
	for (int i = 0; i < size; i++)
	{
		A[i] = new float[size];
		B[i] = new float[size];
		C[i] = new float[size];
		for (int j = 0; j < size; j++)
		{
			A[i][j] = static_cast<float> (rand()) / static_cast<float> (RAND_MAX);
			B[i][j] = static_cast<float> (rand()) / static_cast<float> (RAND_MAX);
			C[i][j] = 0.0f;
		}
	}
	
	//affichage des matrices
	WriteMatrix("A", A, size);
	WriteMatrix("B", B, size);*/
	
	/*
	
	// sequentiel
	float average = 0.0f;
	for (int it = 0; it < num_iter; it++)
	{
		cout << "Algorithme Sequentiel iteration : "<<it << endl;
		t1 = high_resolution_clock::now();

		matrix_multiply(A, B, C, size);

		t2 = high_resolution_clock::now();
		auto duration = std::chrono::duration_cast<std::chrono::microseconds>(t2 - t1).count();
		cout << "Temps d'Execution : " << static_cast <float> (duration) / 1000000.0f << endl;
		average += (static_cast <float> (duration) / 1000000.0f);
	}
	average /= num_iter;
	WriteMatrix("C", C, size);
	WriteData(SEQ_OUTPUT_FILENAME, size, average);

	// parallele
	float parallelAverage;
	for (int n = 2; n < max_threads; n++){
		parallelAverage = 0.0f;
		for (int it = 0; it < num_iter; it++)
		{
			cout << "Algorithme Parallele iteration : " << it << endl;
			t1 = high_resolution_clock::now();

			matrix_multiply_parallel(A, B, C, size, n);

			t2 = high_resolution_clock::now();
			auto duration = std::chrono::duration_cast<std::chrono::microseconds>(t2 - t1).count();
			cout << "Temps d'Execution : " << static_cast <float> (duration) / 1000000.0f << endl;
			parallelAverage += (static_cast <float> (duration) / 1000000.0f);
		}
		parallelAverage /= num_iter;
		WriteData(PARALLEL_OUTPUT_FILENAME, size, parallelAverage, n);
	}
	
		

		float averageBlock = 0.0f;
		for (int it = 0; it < num_iter; it++)
		{
			cout << "Algorithme Sequentiel par block iteration : "<<it << endl;
			t1 = high_resolution_clock::now();

			block_matrix_multiply(A, B, C, size, block_size);

			t2 = high_resolution_clock::now();
			auto duration = std::chrono::duration_cast<std::chrono::microseconds>(t2 - t1).count();
			cout << "Temps d'Execution : " << static_cast <float> (duration) / 1000000.0f << endl;
			averageBlock += (static_cast <float> (duration) / 1000000.0f);
		}
		averageBlock /= num_iter;
		WriteData(SEQ_PER_BLOCK_OUTPUT_FILENAME, size, averageBlock);


		// parallele
		float parallelAverageBlock;
		for (int n = 2; n < max_threads; n++){
			parallelAverageBlock = 0.0f;
			for (int it = 0; it < num_iter; it++)
			{
				cout << "Algorithme Parallele par block iteration : " << it << endl;
				t1 = high_resolution_clock::now();

				block_matrix_multiply_parallel(A, B, C, size, block_size,n);

				t2 = high_resolution_clock::now();
				auto duration = std::chrono::duration_cast<std::chrono::microseconds>(t2 - t1).count();
				cout << "Temps d'Execution : " << static_cast <float> (duration) / 1000000.0f << endl;
				parallelAverageBlock += (static_cast <float> (duration) / 1000000.0f);
			}
			parallelAverageBlock /= num_iter;
			WriteData(PARALLEL_PER_BLOCK_OUTPUT_FILENAME, size, parallelAverageBlock, n);
		}*/

	/*	cout << "Computation has begun(parallel test)" << endl;
		t1 = high_resolution_clock::now();

		test_matrix_mul_parallel(A, B, C, size);

		t2 = high_resolution_clock::now();
		auto duration = std::chrono::duration_cast<std::chrono::microseconds>(t2 - t1).count();
		cout << "Time elapsed: " << static_cast <float> (duration) / 1000000.0f << endl;*/
	
	system("pause");
}



void matrix_multiply(float **A, float **B, float **C, int size)
{
	int i, j, k;
	for (i = 0; i < size; i++)
	{
		for (j = 0; j < size; j++)
		{
			for (k = 0; k < size; k++)
			{
				C[i][j] += A[i][k] * B[k][j];
			}
		}
	}
}

void matrix_multiply_parallel(float **A, float **B, float **C, int size,int nbThreads)
{
	int i, j, k;
	int chunk = 100;
	int tid;
#pragma omp parallel shared(A, B, C, size, chunk) private(k, tid) num_threads(nbThreads)
	{
		tid = omp_get_thread_num();
		if (tid == 0)
		{
			cout << "Number of threads: " << omp_get_num_threads() << endl;
		}
#pragma omp for //schedule (static, chunk)
		for (i = 0; i < size; i++)
		{
			for (j = 0; j < size; j++)
			{
				for (k = 0; k < size; k++)
				{
					C[i][j] += A[i][k] * B[k][j];
				}
			}
		}
	}
}

void block_matrix_multiply(float **A, float **B, float **C, int size, int block_size)
{
	float tmp;
	for (int jj = 0; jj < size; jj += block_size)
	{
		for (int kk = 0; kk < size; kk += block_size)
		{
			for (int i = 0; i < size; i++)
			{
				for (int j = jj; j < ((jj + block_size) > size ? size : (jj + block_size)); j++)
				{
					tmp = 0.0f;
					for (int k = kk; k < ((kk + block_size) > size ? size : (kk + block_size)); k++)
					{
						tmp += A[i][k] * B[k][j];
					}
					C[i][j] += tmp;
				}
			}
		}
	}
}

void block_matrix_multiply_parallel(float **A, float **B, float **C, int size, int block_size,int nbThreads)
{
	int i = 0, j = 0, k = 0, jj = 0, kk = 0;
	float tmp;
	int chunk = 1;
	int tid;
#pragma omp parallel shared(A, B, C, size, chunk) private(i, j, k, jj, kk, tid, tmp) num_threads(nbThreads)
	{
		tid = omp_get_thread_num();
		if (tid == 0)
		{
			cout << "Number of threads: " << omp_get_num_threads() << endl;
		}
#pragma omp for schedule (static, chunk)
		for (jj = 0; jj < size; jj += block_size)
		{
			for (kk = 0; kk < size; kk += block_size)
			{
				for (i = 0; i < size; i++)
				{
					for (j = jj; j < ((jj + block_size) > size ? size : (jj + block_size)); j++)
					{
						tmp = 0.0f;
						for (k = kk; k < ((kk + block_size) > size ? size : (kk + block_size)); k++)
						{
							tmp += A[i][k] * B[k][j];
						}
						C[i][j] += tmp;
					}
				}
			}
		}
	}
}

void display_matrix(float **matrix, int size)
{
	for (int i = 0; i < size; i++)
	{
		for (int j = 0; j < size; j++)
		{
			cout << matrix[i][j] << " ";
		}
		cout << endl;
	}
}

void WriteMatrix(char *fn, float **matrix, int size)
{
	FILE *fp;
	fopen_s(&fp, fn, "w");
	for (int i = 0; i < size; i++)
	{
		for (int j = 0; j < size; j++)
		{
			fprintf(fp, "%f ", matrix[i][j]);
		}
		fprintf(fp,"\n");
	}
	fclose(fp);
}

void WriteData(char *fn, int n, float d)
{
	FILE *fp;
	fopen_s(&fp, fn, "a");
	fprintf(fp, "%d,%lf\n", n, d);
	fclose(fp);
}

void WriteData(char *fn,int n, float d,int ths)
{
	FILE *fp;
	fopen_s(&fp, fn, "a");
	fprintf(fp, "%d,%d,%lf\n", ths,n, d);
	fclose(fp);
}
void CreateData(char *fn)
{
	FILE *fp;
	fopen_s(&fp, fn, "w");
	fprintf(fp, "%s,%s,%s,%s,%s,%s\n", "nbThreads", "size", "executionTime", "speedUp", "efficacity", "realEfficacity");
	fclose(fp);
}
void WriteData(char *fn, int nbThreads, int size, float timeExec,float speedUp,float efficacity ,float realEfficacity)
{
	FILE *fp;
	fopen_s(&fp, fn, "a");
	fprintf(fp, "%d,%d,%lf,%lf,%lf,%lf\n", nbThreads, size, timeExec, speedUp, efficacity, realEfficacity);
	fclose(fp);
}
void scenario(){
	int size, block_size, num_iter, max_threads;
	high_resolution_clock::time_point t1, t2;
	//saisie des parametres
	cin >> size >> block_size >> num_iter >> max_threads;
	//initialisation des matrices
	float **A = new float*[size];
	float **B = new float*[size];
	float **C = new float*[size];

	for (int i = 0; i < size; i++)
	{
		A[i] = new float[size];
		B[i] = new float[size];
		C[i] = new float[size];
		for (int j = 0; j < size; j++)
		{
			A[i][j] = static_cast<float> (rand()) / static_cast<float> (RAND_MAX);
			B[i][j] = static_cast<float> (rand()) / static_cast<float> (RAND_MAX);
			C[i][j] = 0.0f;
		}
	}
	//sequentiel
	float averageBlock = 0.0f;
	for (int it = 0; it < num_iter; it++)
	{
		cout << "Algorithme Sequentiel par block iteration : " << it << endl;
		t1 = high_resolution_clock::now();

		block_matrix_multiply(A, B, C, size, block_size);

		t2 = high_resolution_clock::now();
		auto duration = std::chrono::duration_cast<std::chrono::microseconds>(t2 - t1).count();
		cout << "Temps d'Execution : " << static_cast <float> (duration) / 1000000.0f << endl;
		averageBlock += (static_cast <float> (duration) / 1000000.0f);
	}
	averageBlock /= num_iter;

	// parallele
	float parallelAverageBlock;
	float speedUp;
	float efficacity;
	float realEfficacity;
	for (int n = 2; n < max_threads; n++){
		parallelAverageBlock = 0.0f;
		for (int it = 0; it < num_iter; it++)
		{
			cout << "Algorithme Parallele par block iteration : " << it << endl;
			t1 = high_resolution_clock::now();

			block_matrix_multiply_parallel(A, B, C, size, block_size, n);

			t2 = high_resolution_clock::now();
			auto duration = std::chrono::duration_cast<std::chrono::microseconds>(t2 - t1).count();
			cout << "Temps d'Execution : " << static_cast <float> (duration) / 1000000.0f << endl;
			parallelAverageBlock += (static_cast <float> (duration) / 1000000.0f);
		}
		parallelAverageBlock /= num_iter;
		speedUp = averageBlock / parallelAverageBlock;
		efficacity = speedUp / n;
		realEfficacity = speedUp / 4;
		WriteData(SCENARIO_OUTPUT_FILENAME, n, size, parallelAverageBlock, speedUp, efficacity, realEfficacity);
	}
}