/* 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_before, free_after, total;

    cudaMemGetInfo(&free_before,&total);
    cout << ("%d KB free_before of total %d KB\n",free_before/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);

    cudaMemGetInfo(&free_after,&total);
    cout << ("%d KB free_after of total %d KB\n",free_after/1024,total/1024);

    // Выделено в памяти 2 массива
    int count_of_threads = (free_after - free_before) /  (2 * N * N * sizeof(double));

    cout << "\n" << count_of_threads << endl;

    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;
}