#include <stdio.h>
#include <cublas_v2.h>
#include <thrust/host_vector.h>
#include <thrust/device_vector.h>
#include <thrust/transform.h>
#include <thrust/fill.h>
#include <thrust/sequence.h>
#define N 9000000

struct saxpyFunctor { //Бинарный оператор
    const float a;
    saxpyFunctor(float _a) : a(_a) {}
    __host__ __device__ float operator()(float x, float y) {
        return a * x + y;
    }
};

__global__ void cudaSaxpy(float alpha, float *x, float *y) { //Функция для вычисления Saxpy в CUDA C
    int i = threadIdx.x + blockIdx.x * blockDim.x; //определяем индекс, по координате потока в блоке + координаты блока в сетке * размер блока
    y[i] = alpha * x[i] + y[i]; //вычисляем Saxpy где S - float - тоесть вещественные числа, а axpy это a*x+y (ax plus y = axpy)
}

void saxpy(float a, thrust::device_vector<float>& x, thrust::device_vector<float>& y) { //функция для траст библиотеки
    saxpyFunctor func(a); // инициализация бинарного оператора.
    thrust::transform(x.begin(), x.end(), y.begin(), y.begin(), func); // происходит трансформация, от начала вектора х и У до их конца, соглсно функции func
    //происходит проход от начала x и y, где согласно БИНАРНОЙ функции SAXPY мы выводим результат в операторный вектор y, тоесть мы берем из У информацию, обсчитываем вместе с Х и кладем обратно в ту же ячейку вектора У откуда брали только что данные
}

void CudaKernel(){
    float elapsedTime; // затраченное время
    cudaEvent_t start, stop;  // Объекты событий
    float *x_d, *x_h, *y_h, *y_d; // вещественные указатели на x и y девайса и хоста
    cudaEventCreate(&start); // функция создания события
    cudaEventCreate(&stop);
    cudaMalloc((void**)&x_d, N * sizeof(float)); //выделение глобальной памяти где x_d указатель на память в девайсе, N - размерность
    cudaMalloc((void**)&y_d, N * sizeof(float));
    x_h = (float*)calloc(N, sizeof(float)); //выделение памяти на хосте
    y_h = (float*)calloc(N, sizeof(float));

    for(int i = 0; i < N; i++) { //заполнение векторов на хосте
        x_h[i] = i;
        y_h[i] = 0.87;
    }

    cudaMemcpy(x_d, x_h, N * sizeof(float), cudaMemcpyHostToDevice);    // перенос вектора с хоста на девайс
    cudaMemcpy(y_d, y_h, N * sizeof(float), cudaMemcpyHostToDevice);

    cudaEventRecord(start, 0); // запускаем событие, для начала отсчета таймера
    cudaSaxpy << < N / 256, 256 > >> (3.0, x_d, y_d); //Вызов ядра <<< Число блоков, число нитей в блоке>> ( аргументы функции) туда мы передаем альфу, и векторы на девайсе
    cudaDeviceSynchronize(); // в теории оно тут не нужно
    cudaEventRecord(stop, 0); // фиксируем событие окончания вычислений

    cudaEventSynchronize(stop); //Данная функция ожидает окончание работы всех нитей GPU и прохождение заданного event’а и только потом отдает управление вызывающей программе.
    cudaEventElapsedTime(&elapsedTime, start, stop); //подсчет дельты между событиями
    cudaMemcpy(y_h, y_d, N * sizeof(float), cudaMemcpyDeviceToHost); //копирование с девайса на хост
    printf("CUDA Kernel Time:\n \t %f ms\n", elapsedTime);
    free(x_h); //высвобождение памяти на хосте
    free(y_h);
    cudaFree(y_d); // высвобождение глобальной памяти на девайсе
    cudaFree(x_d);
    cudaEventDestroy(start);
    cudaEventDestroy(stop);
}

void cudaCublas(){
     float elapsedTime;
    cudaEvent_t start, stop;
    float *cx_d, *cx_h, *cy_h, *cy_d;
    cudaEventCreate(&start);
    cudaEventCreate(&stop);

    cudaMallocHost((void**)&cx_h, N * sizeof(float)); //выделение pinned-памяти. особый участок памяти на ОЗУ
    cudaMallocHost((void**)&cy_h, N * sizeof(float));
    cudaMalloc((void**)&cx_d,  N * sizeof(float)); // выделение на девайсе глобальной памяти
    cudaMalloc((void**)&cy_d,  N * sizeof(float));

    for(int i = 0; i < N; i++) { // заполнение на хосте
        cx_h[i] = (float) i;
        cy_h[i] = 0.87f;
    }

    cublasHandle_t cublas_handle; //создание объекта  для работы с библиотекой
    cublasCreate(&cublas_handle); //иницилазация работы с библиотекой

    const int num_rows = N;
    const int stride = 1;

    cublasSetVector(num_rows, sizeof(float), cx_h, stride, cx_d, stride); // инициализация вектора на девайсе
    cublasSetVector(num_rows, sizeof(float), cy_h, stride, cy_d, stride);


    float alpha = 3.0f;

    cudaEventRecord(start, 0);
    cublasSaxpy(cublas_handle, N, &alpha, cx_d, stride, cy_d, stride); // моя функция для высчета суммы векторов
    cudaEventRecord(stop, 0);

    cudaEventSynchronize(stop); // ждем выполнения всех нитей для события завершения
    cudaEventElapsedTime(&elapsedTime, start, stop);

    cublasGetVector(num_rows, sizeof(float), cy_d, stride, cy_h, stride);
    printf("cuBLAS Time:\n \t %f ms\n", elapsedTime);
    cublasDestroy(cublas_handle);
    cudaFreeHost(cx_h);  //высвобождение pinned-памяти
    cudaFreeHost(cy_h);
    cudaFree(cx_d);  // высвобождение хостовой памяти
    cudaFree(cy_d);
    cudaEventDestroy(start);
    cudaEventDestroy(stop);
}


void ThrustLib(){
     float elapsedTime;
    cudaEvent_t start, stop;
    cudaEventCreate(&start);
    cudaEventCreate(&stop);
    thrust::host_vector<float> h1(N); // инициализация векторов на хосте
    thrust::host_vector<float> h2(N);
    thrust::sequence(h1.begin(), h1.end()); // создаем секвенцию в векторе h1 от начала до конца, с заполнением от 1 до N
    thrust::fill(h2.begin(), h2.end(), 0.87); // заполнение вектора неким числом

    thrust::device_vector<float> d1 = h1; // копируем векторы хоста в векторы на девайс
    thrust::device_vector<float> d2 = h2;

    cudaEventRecord(start, 0); //фиксируем событие начала
    saxpy(3.0, d1, d2);
    cudaEventRecord(stop, 0); // кончала

    cudaEventSynchronize(stop); // синхронизируем ждем
    cudaEventElapsedTime(&elapsedTime, start, stop); // подсчитываем финальный итог

    h2 = d2; // обратное копирование
    h1 = d1;

    printf("THRUST Time:\n \t %f ms\n\n", elapsedTime);
    cudaEventDestroy(start);
    cudaEventDestroy(stop);
}


int main() {
    CudaKernel();
    cudaCublas();
    ThrustLib();
    return 0;
}