/* 10.	Вычисление обратной матрицы для заданной матрицы*/


#include <cuda.h>
#include <cuda_runtime.h>
#include "device_launch_parameters.h"
#include <time.h>
#include <iostream>
#include <fstream>
#include <stdlib.h>
#include <stdio.h>
using namespace std;

double *createMass1(double *mass, int x, int y, int N) //создание подматрицы путем вычеркивания столбца и строки
{
	 double *mass1 = new double[N*N];
	int mi = 0;
	for (int i = 0; i < N*N; i++, mi++)
	{
		while (mi / (N + 1) == x || mi % (N + 1) == y)
			mi++;
		mass1[i] = mass[mi];
	}
	return mass1;
}

__device__ double* createMass( double *mass, int x, int y, int N)
{
	double *mass1 = new double[N*N];
	int mi = 0;
	for (int i = 0; i < N*N; i++, mi++)
	{
		while (mi / (N + 1) == x || mi % (N + 1) == y)
			mi++;
		mass1[i] = mass[mi];
	}
	return mass1;
}

double determinant1( double *mass, int N)
{
	 double determ = 1;
	for (int i = 0; i < N*N; i += (N + 1)) //идем по диагонали
	{
		bool d = false;
		for (int j = i; j < N*N; j += i%N)
		{
			 double first = mass[j]; // Выбираем первый слева столбец матрицы
			determ *= first;
			do
			{
				mass[j] /= first; //Все элементы первой строки делят на верхний элемент выбранного столбца
				if (d)
					mass[j] -= mass[i + (j%N - i%N)]; //вычеркивается первая строка и первый столбец
				j++;

			} while (j % N);
			d = true;
		}
	}
	return determ;
}

__device__ double determinant(double *mass, int N)
{
	 double determ = 1;
	for (int i = 0; i < N*N; i += (N + 1)) //идем по диагонали
	{
		bool d = false;
		for (int j = i; j < N*N; j += i%N) // от текущего элемента главной диагонали до последней строки
		{
			 double first = mass[j];// Выбираем первый слева столбец матрицы
			determ *= first;
			do
			{
				mass[j] /= first; //Все элементы первой строки делят на верхний элемент выбранного столбца
				if (d)
					mass[j] -= mass[i + (j%N - i%N)]; //вычеркивается первая строка и первый столбец
				j++;

			} while (j % N);
			d = true;
		}
	}
	return determ;
}

__global__ void division(int i, double m,  double *mass)
{
	mass[threadIdx.x + i] /= m;

}

__global__ void subtraction(int i, int j,  double *mass)
{
	mass[threadIdx.x + i] -= mass[threadIdx.x + j];
}

 double deter(int N,  double *mass)
{
	 double *mass1 = 0;
	// Выделение памяти GPU
	cudaMalloc((void**)&mass1, N * N * sizeof(double));
	// Скопировать входные данные в GPU для обработки
	cudaMemcpy(mass1, mass, N * N * sizeof( double), cudaMemcpyHostToDevice);
	 double determ = 1;
	for (int i = 0; i < N*N; i += (N + 1))
	{
		bool d = false;
		int j, k;
		k = j = i;
		do
		{
			determ *= mass[k]; // Выбираем первый слева столбец матрицы
			division <<<1, (N - i%N) >>>((i + N*(j - i)), mass[k], mass1);
			if (d)
				subtraction <<<1, (N - i%N) >>>((i + N*(j - i)), i, mass1);
			d = true;
			j++;
			k += N;

		} while (j % N);
		// Переписать результаты работы GPU в память CPU:
		cudaMemcpy(mass, mass1, N * N * sizeof( double), cudaMemcpyDeviceToHost);
	}
	return determ;
}

void start1( double *mass, double *temp, int N,  double *d)
{
	 double det = determinant1(d, N);//вычисление детерминанта исходной матрицы
	if (det == 0)
	{
		cout << "Obratnoy matrici ne suchestvuet" << endl;
		return;
	}
	for (int p = 0; p < N*N; p++)
	{
		int i = p;
		int j = i / N;
		i = i % N;
		//создание подматрицы путем вычеркивания i-ого и j-ого столбца для нахождения минора
		 double* A = createMass1(mass, i, j, N - 1);
		//нахождение алгебраического дополнения и деление его на определитель
		temp[p] = (determinant1(A, N - 1) * ((2 + i + j) % 2 ? -1 : 1)) / det;
		delete[]A;
	}
}

__global__ void transp(double *mass, double det, int N)
{
	int k = blockIdx.x * blockDim.x + threadIdx.x;
	mass[k] *= ((2 + (k % N) + (k / N)) % 2 ? -1 : 1);
	mass[k] /= det;
}

__global__ void comp( double *mass, double *temp, int N, double det)
{
	int i = threadIdx.x%N;
	int j = threadIdx.x / N;
	double* A = createMass(mass, i, j, N - 1);
	temp[threadIdx.x] = determinant(A, N - 1);
	delete[]A;
}


void start( double *mass, double *temp, int N, double *d)
{
	size_t free, total;
	printf("\n");

	cudaMemGetInfo(&free,&total);

	printf("%d KB free of total %d KB\n",free/1024,total/1024);

	if (N >= 32)
    {
        cout << "Matrica vichodit za granici pamyati GPU adaptera" << endl;
        return;
    }

	double det = deter(N, d);
	if (det == 0)
	{
		cout << "Obratnoy matrici ne suchestvuet" << endl;
		return;
	}
	double *temp1 = 0, *mass1 = 0;
	// Выделение памяти GPU
	cudaMalloc((void**)&temp1, N * N * sizeof(double));
	cudaMalloc((void**)&mass1, N * N * sizeof(double));
	// Скопировать входные данные в GPU для обработки
	cudaMemcpy(mass1, mass, N * N * sizeof(double), cudaMemcpyHostToDevice);
	// Инициализация ядра в GPU
	comp<<< 1, N*N >>>(mass1, temp1, N, det);
	transp<<< 1, N*N >>>(temp1, det, N);
	// Переписать результаты работы GPU в память CPU:
	cudaMemcpy(temp, temp1, N * N * sizeof(double), cudaMemcpyDeviceToHost);
	cudaFree(temp1);
	cudaFree(mass1);
}

bool compare(double *mass1, double *mass2, int N)
{
	for (int i = 0; i < N; i++)
	{
		if (mass1[i] != mass2[i])
			return false;
	}
	return true;
}

int main()
{
	setlocale(LC_CTYPE, "rus");
	srand(time(0));
	int N;
	cout << "Vvedite razmer matrici: ";
	cin >> N;
	double *mass = new double[N*N];
	double *temp = new double[N*N];
	double *_mass = new double[N*N];
	double *_temp = new double[N*N];
	double *d = new double[N*N];
	double *_d = new double[N*N];

	for (int i = 0; i < N*N; i++)
		mass[i] = _mass[i] = d[i] = _d[i] = ((double)rand()) / 100000;



	double t1 = clock();
	start(mass, temp, N, d);
	t1 = clock() - t1;
	cout << endl << endl << "Parallel: " << t1 / CLOCKS_PER_SEC << " sec." << endl;

	double t = clock();
	start1(_mass, _temp, N, _d);
	t = clock() - t;
	cout << endl << "Posledovat: " << t / CLOCKS_PER_SEC << " sec." << endl;

	if (t - t1 > 0)
		cout << endl << endl << "Parallel bystree posledovat v " << int(t / t1) << " raz." << endl;
	else
		cout << endl << endl << "Parallel medlennee posledovat v " << int(t1 / t) << " raz." << endl;

	if (compare(temp, _temp, N))
		cout << endl << "Resultat identichni!" << endl << endl;
	else
		cout << endl << "Resultat razlichni!" << endl << endl;

	delete[]temp;
	delete[]mass;
	delete[]_temp;
	delete[]_mass;
}