Шпора на ПГП

#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

//SAXPY — это комбинация скалярного умножения и сложения векторов, то есть: в качестве входных данных
//принимается два 32-разрядных вектора X и Y с N элементами каждый и скалярное значение A. Далее происходит
//умножение каждого элемента X [i] на A, а получение результата приходит к Y [i].

//Функтор - объект, который можно использовать как функцию
//Функтор получающий на вход два аргумента и возвращающий результат их обработки как переопределенный оператор ().
struct saxpyFunctor {//Бинарный оператор
    const float a;
    saxpyFunctor(float _a) : a(_a) {}
    __host__ __device__ float operator()(float x, float y) {
        return a * x + y;
    }
};

//threadIdx.x это потоки в блоке
//blockIdx.x это сам блок
//blockDim.x это размерность блока
//представим что у нас есть цикл и потоки допустим на блок у нас от 1 до 5
//блоки исчисляются с 0 и далее а размерность блока допустим 5
//тогда первый цикл первый поток + 0 блок * 5 = 1 индекс в виюдюхе
//второй поток + 0 блок *5 = 2 индекс
//потом допустим блок мы этот забили таким образом, значит переменная blockIdx.x + 1 сделала
//получаем 1 поток в блоке + 1 блок *5 = 6 индекс в видюхе
__global__ void cudaSaxpy(float alpha, float* x, float* y) {
    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);
    //выделение глобальной памяти где x_d указатель на память в девайсе, N - размерность
    cudaMalloc((void**)&x_d, N * sizeof(float));
    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);// фиксируем событие окончания вычислений
    //Данная функция ожидает окончание работы всех нитей GPU и прохождение заданного event’а и только потом отдает управление вызывающей программе.
    cudaEventSynchronize(stop);//Ожидание до завершения всех работ, захваченных в данный момент в события
    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);
    //выделение закрепленной (pinned) памяти, особый участок памяти на ОЗУ
    cudaMallocHost((void**)&cx_h, N * sizeof(float));
    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. Контекст библиотеки cuBLAS должен быть
    //инициализирован с использованием cublasCreate () и возвращенный дескриптор должен быть передан всем последующим вызовам
    //библиотечных функций. Контекст должен быть уничтожен в конце, используяcublasDestroy (),
    cublasHandle_t cublas_handle;//создание объекта  для работы с библиотекой
    cublasCreate(&cublas_handle);//иницилазация работы с библиотекой

    const int num_rows = N;
    const int num_cols = 1;
    const int stride = 1;
    //const size_t elem_size = sizeof(float);

    //Эта функция копирует плитку строки х столбцы элементы из матрицы A в пространстве памяти хоста в матрицу Вв пространстве памяти GPU
    //cublasSetMatrix(num_rows, num_cols, elem_size, cx_h, num_rows, cx_d, num_rows);
    //cublasSetMatrix(num_rows, num_cols, elem_size, cy_h, num_rows, cy_d, num_rows);

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

    const int stride = 1;
    float alpha = 3.0f;

    cudaEventRecord(start, 0);
    //Сложение векторов (реальная одинарная точность).
    //cublas_handle - контекст CUBLAS.
    //N - количество элементов во входных векторах.
    //alpha ( numpy.float32 ) - скаляр .
    //cx_d - указатель на входной вектор с одинарной точностью.
    //stride -  Интервал хранения между элементами x .
    //cy_d - указатель на вектор ввода / вывода одинарной точности.
    //stride - Интервал хранения между элементами y .
    cublasSaxpy(cublas_handle, N, &alpha, cx_d, stride, cy_d, stride);
    cudaEventRecord(stop, 0);

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

    //Эта функция копирует плитку строки х столбцы элементы из матрицы A в пространстве памяти GPU до матрицы Вв пространстве памяти хоста
    //cublasGetMatrix(num_rows, num_cols, elem_size, cx_d, num_rows, cx_h, num_rows);
    //cublasGetMatrix(num_rows, num_cols, elem_size, cy_d, num_rows, cy_h, num_rows);

    cublasGetVector(num_rows, sizeof(float), cy_d, stride, cy_h, stride);
    printf("cuBLAS Time:\n \t %f ms\n", elapsedTime);

    //Эта функция высвобождает аппаратные ресурсы, используемые библиотекой cuBLAS
    cublasDestroy(cublas_handle);

    //Освобождает пространство памяти, на которое указывает hostPtr, которое должно быть возвращено предыдущим вызовом cudaMallocHost () или cudaHostAlloc ()
    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);

    //контейнеры, параметризованные типом хранящихся в них значений
    //хранит содержимое вектора в памяти CPU
    thrust::host_vector<float> h1(N);// инициализация векторов на хосте
    thrust::host_vector<float> h2(N);

    //заполняет диапазон [first, last)последовательностью чисел.
    //Для каждого итератора iв диапазоне [first, last)эта версия sequenceвыполняет присваивание *i = init + (i - first).
    thrust::sequence(h1.begin(), h1.end());// создаем секвенцию в векторе h1 от начала до конца, с заполнением от 1 до N

    //fillприсваивает значение value каждому элементу в диапазоне [first, last)
    thrust::fill(h2.begin(), h2.end(), 0.87);// заполнение вектора неким числом

    //контейнеры, параметризованные типом хранящихся в них значений
    //хранит содержимое вектора в памяти GPU
    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;
}