Еще одна статья о CUDA — зачем?
На Хабре было уже немало хороших статей по CUDA — раз, два и другие. Однако, поиск комбинации «CUDA scan» выдал всего 2 статьи никак не связанные с, собственно, алгоритмом scan на GPU — а это один из самых базовых алгоритмов. Поэтому, вдохновившись только что просмотренным курсом на Udacity — Intro to Parallel Programming, я и решился написать более полную серию статей о CUDA. Сразу скажу, что серия будет основываться именно на этом курсе, и если у вас есть время — намного полезнее будет пройти его.
Содержание
На данный момент планируются следующие статьи:
Часть 1: Введение.
Часть 2: Аппаратное обеспечение GPU и шаблоны параллельной коммуникации.
Часть 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 вы просто пишете код на своем любимом языке программирования (список поддерживаемых языков, не учитывая С и С++), после чего компилятор CUDA сгенерирует код отдельно для хоста и отдельно для устройства. Небольшая оговорка: код для устройства должен быть написан только на языке C с некоторыми 'CUDA-расширениями'.
Основные этапы CUDA-программы
- Хост выделяет нужное количество памяти на устройстве.
- Хост копирует данные из своей памяти в память устройства.
- Хост стартует выполнение определенных ядер на устройстве.
- Устройство выполняет ядра.
- Хост копирует результаты из памяти устройства в свою память.
Естественно, для наибольшей эффективности использования GPU нужно чтобы соотношение времени, потраченного на работу ядер, к времени, потраченного на выделение памяти и перемещение данных, было как можно больше.
Ядра
Рассмотрим более детально процесс написания кода для ядер и их запуска. Важный принцип — ядра пишутся как (практически) обычные последовательные программы — то-есть вы не увидите создания и запуска потоков в коде самих ядер. Вместо этого, для организации параллельных вычислений GPU запустит большое количество копий одного и того же ядра в разных потоках — а точнее, вы сами говорите сколько потоков запустить. И да, возвращаясь к вопросу эффективности использования GPU — чем больше потоков вы запускаете (при условии что все они будут выполнять полезную работу) — тем лучше.
Код для ядер отличается от обычного последовательного кода в таких моментах:
- Внутри ядер вы имеете возможность узнать «идентификатор» или, проще говоря, позицию потока, который сейчас выполняется — используя эту позицию мы добиваемся того, что одно и то же ядро будет работать с разными данными в зависимости от потока, в котором оно запущено. Кстати, такая организация параллельных вычислений называется SIMD (Single Instruction Multiple Data) — когда несколько процессоров выполняют одновременно одну и ту же операцию но на разных данных.
- В некоторых случаях в коде ядра необходимо использовать различные способы синхронизации.
Каким же образом мы задаем количество потоков, в которых будет запущено ядро? Поскольку 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. Также, для простоты работы с файлами изображений будем использовать OpenCV, а для сравнения производительности CPU и GPU — OpenMP.
Задачу поставим довольно простую: конвертация цветного изображения в оттенки серого. Для этого, яркость пиксела 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 к однопоточному коду — в этом вся красота и мощь OpenMP. Я пробовал поиграться с параметром schedule, но получалось только хуже, чем без него.
Наконец, переходим к 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-репозитория 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.
Автор: VladGorbatiuk
