Параллельное программирование с CUDA. Часть 1: Введение

Еще одна статья о CUDA — зачем?

На Хабре было уже немало хороших статей по CUDA — раз ^[1], два ^[2] и другие. Однако, поиск комбинации «CUDA scan ^[3]» выдал всего 2 статьи никак не связанные с, собственно, алгоритмом scan на GPU — а это один из самых базовых алгоритмов. Поэтому, вдохновившись только что просмотренным курсом на Udacity — Intro to Parallel Programming ^[4], я и решился написать более полную серию статей о CUDA. Сразу скажу, что серия будет основываться именно на этом курсе, и если у вас есть время — намного полезнее будет пройти его.

Содержание

На данный момент планируются следующие статьи:
Часть 1: Введение.
Часть 2: Аппаратное обеспечение GPU и шаблоны параллельной коммуникации. ^[5]
Часть 3: Фундаментальные алгоритмы GPU: свертка (reduce), сканирование (scan) и гистограмма (histogram).
Часть 4: Фундаментальные алгоритмы GPU: уплотнение (compact), сегментированное сканирование (segmented scan), сортировка. Практическое применение некоторых алгоритмов.
Часть 5: Оптимизация GPU программ.
Часть 6: Примеры параллелизации последовательных алгоритмов.
Часть 7: Дополнительные темы параллельного программирования, динамический параллелизм.

Задержка vs пропускная способность

Первый вопрос, который должен задать каждый перед применением GPU для решения своих задач — а для каких целей хорош GPU, когда стоит его применять? Для ответа нужно определить 2 понятия:
Задержка (latency) — время, затрачиваемое на выполнение одной инструкции/операции.
Пропускная способность — количество инструкций/операций, выполняемых за единицу времени.
Простой пример: имеем легковой автомобиль со скоростью 90 км/ч и вместимостью 4 человека, и автобус со скоростью 60 км/ч и вместимостью 20 человек. Если за операцию принять перемещение 1 человека на 1 километр, то задержка легкового автомобиля — 3600/90=40с — за столько секунд 1 человек преодолеет расстояние в 1 километр, пропускная способность автомобиля — 4/40=0.1 операций/секунду; задержка автобуса — 3600/60=60с, пропускная способность автобуса — 20/60=0.3(3) операций/секунду.
Так вот, CPU — это автомобиль, GPU — автобус: он имеет большую задержку но также и большую пропускную способность. Если для вашей задачи задержка каждой конкретной операции не настолько важна как количество этих операций в секунду — стоит рассмотреть применение GPU.

Базовые понятия и термины CUDA

Итак, разберемся с терминологией CUDA:

• Устройство (device) — GPU. Выполняет роль «подчиненного» — делает только то, что ему говорит CPU.
• Хост (host) — CPU. Выполняет управляющую роль — запускает задачи на устройстве, выделяет память на устройстве, перемещает память на/с устройства. И да, использование CUDA предполагает, что как устройство так и хост имеют свою отдельную память.
• Ядро (kernel) — задача, запускаемая хостом на устройстве.

При использовании CUDA вы просто пишете код на своем любимом языке программирования (список ^[6] поддерживаемых языков, не учитывая С и С++), после чего компилятор CUDA сгенерирует код отдельно для хоста и отдельно для устройства. Небольшая оговорка: код для устройства должен быть написан только на языке C с некоторыми 'CUDA-расширениями'.

Основные этапы CUDA-программы

1. Хост выделяет нужное количество памяти на устройстве.
2. Хост копирует данные из своей памяти в память устройства.
3. Хост стартует выполнение определенных ядер на устройстве.
4. Устройство выполняет ядра.
5. Хост копирует результаты из памяти устройства в свою память.

Естественно, для наибольшей эффективности использования GPU нужно чтобы соотношение времени, потраченного на работу ядер, к времени, потраченного на выделение памяти и перемещение данных, было как можно больше.

Ядра

Рассмотрим более детально процесс написания кода для ядер и их запуска. Важный принцип — ядра пишутся как (практически) обычные последовательные программы — то-есть вы не увидите создания и запуска потоков в коде самих ядер. Вместо этого, для организации параллельных вычислений GPU запустит большое количество копий одного и того же ядра в разных потоках — а точнее, вы сами говорите сколько потоков запустить. И да, возвращаясь к вопросу эффективности использования GPU — чем больше потоков вы запускаете (при условии что все они будут выполнять полезную работу) — тем лучше.
Код для ядер отличается от обычного последовательного кода в таких моментах:

1. Внутри ядер вы имеете возможность узнать «идентификатор» или, проще говоря, позицию потока, который сейчас выполняется — используя эту позицию мы добиваемся того, что одно и то же ядро будет работать с разными данными в зависимости от потока, в котором оно запущено. Кстати, такая организация параллельных вычислений называется SIMD ^[7] (Single Instruction Multiple Data) — когда несколько процессоров выполняют одновременно одну и ту же операцию но на разных данных.
2. В некоторых случаях в коде ядра необходимо использовать различные способы синхронизации.

Каким же образом мы задаем количество потоков, в которых будет запущено ядро? Поскольку GPU это все таки Graphics Processing Unit, то это, естественно, повлияло на модель CUDA, а именно на способ задания количества потоков:

• Сначала задаются размеры так называемой сетки (grid), в 3D координатах: grid_x, grid_y, grid_z. В результате, сетка будет состоять из grid_x*grid_y*grid_z блоков.
• Потом задаются размеры блока в 3D координатах: block_x, block_y, block_z. В результате, блок будет состоять из block_x*block_y*block_z потоков. Итого, имеем grid_x*grid_y*grid_z*block_x*block_y*block_z потоков. Важное замечание — максимальное количество потоков в одном блоке ограничено и зависит от модели GPU — типичны значения 512 (более старые модели) и 1024 (более новые модели).
• Внутри ядра доступны переменные threadIdx и blockIdx с полями x, y, z — они содержат 3D координаты потока в блоке и блока в сетке соответственно. Также доступны переменные blockDim и gridDim с теми же полями — размеры блока и сетки соответственно.

Как видите, данный способ запуска потоков действительно подходит для обработки 2D и 3D изображений: например, если нужно определенным образом обработать каждый пиксел 2D либо 3D изображения, то после выбора размеров блока (в зависимости от размеров картинки, способа обработки и модели GPU) размеры сетки выбираются такими, чтобы было покрыто все изображение, возможно, с избытком — если размеры изображения не делятся нацело на размеры блока.

Пишем первую программу на CUDA

Довольно теории, время писать код. Инструкции по установке и конфигурации CUDA для разных ОС — docs.nvidia.com/cuda/index.html ^[8]. Также, для простоты работы с файлами изображений будем использовать OpenCV ^[9], а для сравнения производительности CPU и GPU — OpenMP ^[10].
Задачу поставим довольно простую: конвертация цветного изображения в оттенки серого ^[11]. Для этого, яркость пиксела pix в серой шкале считается по формуле: Y = 0.299*pix.R + 0.587*pix.G + 0.114*pix.B.
Сначала напишем скелет программы:

#include <chrono>
#include <iostream>
#include <cstring>
#include <string>

#include <opencv2/core/core.hpp>
#include <opencv2/highgui/highgui.hpp>
#include <opencv2/opencv.hpp>
#include <vector_types.h>

#include "openMP.hpp"
#include "CUDA_wrappers.hpp"
#include "common/image_helpers.hpp"

using namespace cv;
using namespace std;

int main( int argc, char** argv )
{
  using namespace std::chrono;

  if( argc != 2)
  {
    cout <<" Usage: convert_to_grayscale imagefile" << endl;
    return -1;
  }

  Mat image, imageGray;
  uchar4 *imageArray;
  unsigned char *imageGrayArray;

  prepareImagePointers(argv[1], image, &imageArray, imageGray, &imageGrayArray, CV_8UC1);

  int numRows = image.rows, numCols = image.cols;

  auto start = system_clock::now();
  RGBtoGrayscaleOpenMP(imageArray, imageGrayArray, numRows, numCols);
  auto duration = duration_cast<milliseconds>(system_clock::now() - start);
  cout<<"OpenMP time (ms):" << duration.count() << endl;

  memset(imageGrayArray, 0, sizeof(unsigned char)*numRows*numCols);  

  RGBtoGrayscaleCUDA(imageArray, imageGrayArray, numRows, numCols);

  return 0;
}

Тут все довольно очевидно — читаем файл с изображением, подготавливаем указатели на цветное и в оттенках серого изображение, запускаем вариант
с OpenMP и вариант с CUDA, замеряем время. Функция prepareImagePointers имеет следующий вид:

template <class T1, class T2>
void prepareImagePointers(const char * const inputImageFileName,
                          cv::Mat& inputImage, 
                          T1** inputImageArray, 
                          cv::Mat& outputImage,
                          T2** outputImageArray, 
                          const int outputImageType)
{
  using namespace std;
  using namespace cv;

  inputImage = imread(inputImageFileName, IMREAD_COLOR);

  if (inputImage.empty()) 
  {
    cerr << "Couldn't open input file." << endl;
    exit(1);
  }

  //allocate memory for the output
  outputImage.create(inputImage.rows, inputImage.cols, outputImageType);

  cvtColor(inputImage, inputImage, cv::COLOR_BGR2BGRA);

  *inputImageArray = (T1*)inputImage.ptr<char>(0);
  *outputImageArray  = (T2*)outputImage.ptr<char>(0); 
}

Я пошел на небольшую хитрость: дело в том, что мы выполняем очень мало работы на каждый пиксел изображения — то-есть при варианте с CUDA встает упомянутая выше проблема соотношения времени выполнения полезных операций к времени выделения памяти и копирования данных, и в результате общее время CUDA варианта будет больше OpenMP варианта, а мы же хотим показать что CUDA быстрее:) Поэтому для CUDA будет измеряться только время, потраченное на выполнение собственно конвертации изображения — без учета операций с памятью. В свое оправдание скажу, что для большого класса задач время полезной работы будет все-таки доминировать, и CUDA будет быстрее даже с учетом операций с памятью.
Далее напишем код для OpenMP варианта:

#include <stdio.h>

#include <omp.h>
#include <vector_types.h>

void RGBtoGrayscaleOpenMP(uchar4 *imageArray, unsigned char *imageGrayArray, int numRows, int numCols)
{
    #pragma omp parallel for collapse(2)
    for (int i = 0; i < numRows; ++i)
    {
        for (int j = 0; j < numCols; ++j)
        {
            const uchar4 pixel = imageArray[i*numCols+j];
            imageGrayArray[i*numCols+j] = 0.299f*pixel.x + 0.587f*pixel.y+0.114f*pixel.z;
        }
    }
}

Все довольно прямолинейно — мы всего лишь добавили директиву omp parallel for ^[12] к однопоточному коду — в этом вся красота и мощь OpenMP. Я пробовал поиграться с параметром schedule ^[13], но получалось только хуже, чем без него.
Наконец, переходим к CUDA. Тут распишем более детально. Сначала нужно выделить память под входные данные, переместить их с CPU на GPU и выделить память под выходные данные:

void RGBtoGrayscaleCUDA(const uchar4 * const h_imageRGBA, unsigned char* const h_imageGray, size_t numRows, size_t numCols)
{
  uchar4 *d_imageRGBA;
  unsigned char *d_imageGray;
  const size_t numPixels = numRows * numCols;
  cudaSetDevice(0);
  checkCudaErrors(cudaGetLastError());
  //allocate memory on the device for both input and output
  checkCudaErrors(cudaMalloc(&d_imageRGBA, sizeof(uchar4) * numPixels));
  checkCudaErrors(cudaMalloc(&d_imageGray, sizeof(unsigned char) * numPixels));

  //copy input array to the GPU
  checkCudaErrors(cudaMemcpy(d_imageRGBA, h_imageRGBA, sizeof(uchar4) * numPixels, cudaMemcpyHostToDevice));

Стоит обратить внимание на стандарт именования переменных в CUDA — данные на CPU начинаются с h_ (host), данные да GPU — с d_ (device). checkCudaErrors — макрос, взят с github-репозитория ^[14] Udacity курса. Имеет следующий вид:

#include <cuda.h>

#define checkCudaErrors(val) check( (val), #val, __FILE__, __LINE__)

template<typename T>
void check(T err, const char* const func, const char* const file, const int line) {
  if (err != cudaSuccess) {
    std::cerr << "CUDA error at: " << file << ":" << line << std::endl;
    std::cerr << cudaGetErrorString(err) << " " << func << std::endl;
    exit(1);
  }
}

cudaMalloc — аналог malloc для GPU, cudaMemcpy — аналог memcpy, имеет дополнительный параметр в виде enum-а, который указывает тип копирования: cudaMemcpyHostToDevice, cudaMemcpyDeviceToHost, cudaMemcpyDeviceToDevice.
Далее необходимо задать размеры сетки и блока и вызвать ядро, не забыв измерить время:

 dim3 blockSize;
  dim3 gridSize;
  int threadNum;

  cudaEvent_t start, stop;
  cudaEventCreate(&start);
  cudaEventCreate(&stop);

  threadNum = 1024;
  blockSize = dim3(threadNum, 1, 1);
  gridSize = dim3(numCols/threadNum+1, numRows, 1);
  cudaEventRecord(start);
  rgba_to_grayscale_simple<<<gridSize, blockSize>>>(d_imageRGBA, d_imageGray, numRows, numCols);
  cudaEventRecord(stop);
  cudaEventSynchronize(stop);
  cudaDeviceSynchronize(); checkCudaErrors(cudaGetLastError());
  float milliseconds = 0;
  cudaEventElapsedTime(&milliseconds, start, stop);
  std::cout << "CUDA time simple (ms): " << milliseconds << std::endl;

Обратите внимание на формат вызова ядра — kernel_name<<<gridSize, blockSize>>>. Код самого ядра также не очень сложный:

__global__
void rgba_to_grayscale_simple(const uchar4* const d_imageRGBA,
                              unsigned char* const d_imageGray,
                              int numRows, int numCols)
{
    int y = blockDim.y*blockIdx.y + threadIdx.y;
    int x = blockDim.x*blockIdx.x + threadIdx.x;
    if (x>=numCols || y>=numRows)
      return;
    const int offset = y*numCols+x;
    const uchar4 pixel = d_imageRGBA[offset];
    d_imageGray[offset] = 0.299f*pixel.x + 0.587f*pixel.y+0.114f*pixel.z;
}

Здесь мы вычисляем координаты y и x обрабатываемого пиксела, используя ранее описанные переменные threadIdx, blockIdx и blockDim, ну и выполняем конвертацию. Обратите внимание на проверку if (x>=numCols || y>=numRows) — так как размеры изображения не обязательно будут делится нацело на размеры блоков, некоторые блоки могут «выходить за рамки» изображения — поэтому необходима эта проверка. Также, функция ядра должна помечаться спецификатором __global__ .
Последний шаг — cкопировать результат назад с GPU на CPU и освободить выделенную память:

  checkCudaErrors(cudaMemcpy(h_imageGray, d_imageGray, sizeof(unsigned char) * numPixels, cudaMemcpyDeviceToHost));
  cudaFree(d_imageGray);
  cudaFree(d_imageRGBA);

Кстати, CUDA позволяет использовать C++ компилятор для host-кода — так что запросто можно написать обертки для автоматического освобождения памяти.
Итак, запускаем, измеряем:

OpenMP time (ms):45
CUDA time simple (ms): 43.1941

Получилось как-то не очень впечатляюще:) А проблема все та же — слишком мало работы выполняется над каждым пикселом — мы запускаем тысячи потоков, каждый из которых отрабатывает практически моментально. В случае с CPU такой проблемы не возникает — OpenMP запустит сравнительно малое количество потоков (8 в моем случае) и разделит работу между ними поровну — таким образом процессоры будет занят практически на все 100%, в то время как с GPU мы, по сути, не используем всю его мощь. Решение довольно очевидное — обрабатывать несколько пикселов в ядре. Новое, оптимизированное, ядро будет выглядеть следующим образом:

#define WARP_SIZE 32

__global__
void rgba_to_grayscale_optimized(const uchar4* const d_imageRGBA,
                                 unsigned char* const d_imageGray,
                                 int numRows, int numCols,
                                 int elemsPerThread)
{
    int y = blockDim.y*blockIdx.y + threadIdx.y;
    int x = blockDim.x*blockIdx.x + threadIdx.x;
    const int loop_start =  (x/WARP_SIZE * WARP_SIZE)*(elemsPerThread-1)+x;
    for (int i=loop_start, j=0; j<elemsPerThread && i<numCols; i+=WARP_SIZE, ++j)
    {
      const int offset = y*numCols+i;
      const uchar4 pixel = d_imageRGBA[offset];
      d_imageGray[offset] = 0.299f*pixel.x + 0.587f*pixel.y+0.114f*pixel.z;
    }
}

Здесь не все так просто как с предыдущим ядром. Если разобраться, теперь каждый поток будет обрабатывать elemsPerThread пикселов, причем не подряд, а с расстоянием в WARP_SIZE между ними. Что такое WARP_SIZE, почему оно равно 32, и зачем обрабатывать пиксели пободным образом, будет более детально рассказано в следующих частях, сейчас только скажу что этим мы добиваемся более эффективной работы с памятью. Каждый поток теперь обрабатывает elemsPerThread пикселов с расстоянием в WARP_SIZE между ними, поэтому x-координата первого пиксела для этого потока исходя из его позиции в блоке теперь рассчитывается по несколько более сложной формуле чем раньше.
Запускается это ядро следующим образом:

  threadNum=128;
  const int elemsPerThread = 16;
  blockSize = dim3(threadNum, 1, 1);
  gridSize = dim3(numCols / (threadNum*elemsPerThread) + 1, numRows, 1);
  cudaEventRecord(start);
  rgba_to_grayscale_optimized<<<gridSize, blockSize>>>(d_imageRGBA, d_imageGray, numRows, numCols, elemsPerThread);
  cudaEventRecord(stop);
  cudaEventSynchronize(stop);
  cudaDeviceSynchronize(); checkCudaErrors(cudaGetLastError());
  milliseconds = 0;
  cudaEventElapsedTime(&milliseconds, start, stop);
  std::cout << "CUDA time optimized (ms): " << milliseconds << std::endl;

Количество блоков по x-координате теперь рассчитывается как numCols / (threadNum*elemsPerThread) + 1 вместо numCols / threadNum + 1. В остальном все осталось так же.
Запускаем:

OpenMP time (ms):44
CUDA time simple (ms): 53.1625
CUDA time optimized (ms): 15.9273

Получили прирост по скорости в 2.76 раза (опять же, не учитывая время на операции с памятью) — для такой простой проблемы это довольно неплохо. Да-да, эта задача слишком простая — с ней достаточно хорошо справляется и CPU. Как видно из второго теста, простая реализация на GPU может даже проигрывать по скорости реализации на CPU.
На сегодня все, в следующей части рассмотрим аппаратное обеспечение GPU и основные шаблоны параллельной коммуникации.
Весь исходный код доступен на bitbucket ^[15].

Автор: VladGorbatiuk

Источник ^[16]