- PVSM.RU - https://www.pvsm.ru -
Часть 1: Введение. [1]
Часть 2: Аппаратное обеспечение GPU и шаблоны параллельной коммуникации.
Часть 3: Фундаментальные алгоритмы GPU: свертка (reduce), сканирование (scan) и гистограмма (histogram).
Часть 4: Фундаментальные алгоритмы GPU: уплотнение (compact), сегментированное сканирование (segmented scan), сортировка. Практическое применение некоторых алгоритмов.
Часть 5: Оптимизация GPU программ.
Часть 6: Примеры параллелизации последовательных алгоритмов.
Часть 7: Дополнительные темы параллельного программирования, динамический параллелизм.
Что собой представляют параллельные вычисления? Не что иное, как множество потоков, решающие определенную задачу кооперативно. Ключевое слово здесь «кооперативно» — для достижения кооперации нужно применять определенные механизмы коммуникации между потоками. При использовании CUDA, коммуникация осуществляется через память: потоки могут читать входные данные, изменять выходные данные либо обмениваться «промежуточными» результатами.
В зависимости от того, каким именно образом потоки общаются через память, выделяют различные шаблоны параллельной коммуникации.
В предыдущей части в качестве простого примера использования CUDA была рассмотрена задача конвертации цветного изображения в оттенки серого. Для этого, интенсивность каждого пиксела выходного изображения в оттенках серого рассчитывалась по формуле I=A*pix.R+B*pix.G+C*pix.B, где A, B, C — константы, pix — соответствующий пиксел исходного изображения. Графически этот процесс выглядит следующим образом:
Если абстрагироваться от самого способа расчета выходного значения, получаем первый шаблон параллельной коммуникации — map: над каждым элементом входных данных, с индексом i, выполняется одна и та же функция, а результат сохраняется в выходном массиве данных под тем же индексом i. Шаблон map очень эффективен на GPU, и к тому же просто выражается в рамках CUDA — достаточно запустить по одному потоку на каждый входной элемент (как и было сделано в предыдущей части для задачи конвертации изображения). Однако, лишь малую часть задач можно решить используя только этот шаблон.
Если же для вычисления выходного значения с индексом i используются несколько входных элементов, то такой шаблон называют gather, и выглядеть он может так:
Или так:
Эффективность реализации данного шаблона на CUDA зависит от того, какие именно входные значения используются при расчете выходного и их количества — лучше всего когда используется небольшое количество идущих подряд элементов.
Обратный шаблон, scatter — каждый входной элемент влияет на несколько (либо один) выходных элементов, графически выглядит так же как и gather, однако меняется смысл: теперь мы «отталкиваемся» не от выходных элементов, для которых рассчитываем их значение, а от входных, которые влияют на значения определенных выходных элементов. Довольно часто одна и та же задача может быть решена в рамках как шаблона gather, так и scatter. Например, если мы хотим усреднить 3 соседних входных элемента и записать их в выходной массив, то можно:
Как вы скорее всего догадались, при использовании подхода scatter встает проблема синхронизации, так как несколько потоков могут пытаться модифицировать один и тот же выходной элемент одновременно.
Также стоит выделить подвид шаблона gather — stencil: в данном шаблоне накладывается ограничение на входные элементы, которые участвуют в вычислении выходного элемента — а именно, это могут быть только соседние элементы. В случае с 2D/3D изображениями, могут использоваться различные виды данного шаблона, например двумерный трафарет фон Неймана:
или двумерный трафарет Мура:
В связи с этим ограничением, шаблон stencil обычно довольно эффективно реализуется в рамках CUDA: достаточно запустить по одному потоку на каждый выходной элемент, при этом поток будет считывать нужные ему входные элементы. При такой организации вычислений, эффективность обеспечивается двумя факторами:
Рассмотрим высокоуровневую структуру аппаратного обеспечения GPU:
Согласно модели CUDA, программист разбивает задачу на блоки, а блоки на потоки. Каким же образом выполняется сопоставление этих программных сущностей с выше описанными аппаратными блоками GPU?
Исходя из всего этого, CUDA предоставляет следующие гарантии:
CUDA не гарантирует что:
Итак, перечислим основные механизмы синхронизации, предоставляемые CUDA:
if (threadIdx.x % 2 == 0)
{
...
}
...
, то половина потоков варпа (те, которые с нечетным индексом) будут ждать, пока вторая половина выполнит код, находящийся внутри if-а. Поэтому следует избегать таких ситуаций.
Переходим к практике. В качестве примера, иллюстрирующего изложенную теорию, напишем программу, выполняющую размытие изображения по Гауссу [4]. Принцип работы следующий: значение каналов R, G, B пиксела в выходном размытом изображении рассчитывается как взвешенная сумма значений каналов R, G, B (соответственно) всех пикселов исходного изображения в некотором трафарете:
Значения весов рассчитываются используя 2D распределение Гаусса, но как именно это делается не слишком важно для нашей задачи.
Как можно увидеть из описания задачи, для реализации данного алгоритма вполне естественно выбрать шаблон stencil — ведь каждый пиксел выходного изображения рассчитывается исходя из соответствующих соседних пикселей исходного изображения.
Начнем со скелета программы:
#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"
void prepareFilter(float **filter, int *filterWidth, float *filterSigma)
{
static const int blurFilterWidth = 9;
static const float blurFilterSigma = 2.;
*filter = new float[blurFilterWidth * blurFilterWidth];
*filterWidth = blurFilterWidth;
*filterSigma = blurFilterSigma;
float filterSum = 0.f;
const int halfWidth = blurFilterWidth/2;
for (int r = -halfWidth; r <= halfWidth; ++r)
{
for (int c = -halfWidth; c <= halfWidth; ++c)
{
float filterValue = expf( -(float)(c * c + r * r) / (2.f * blurFilterSigma * blurFilterSigma));
(*filter)[(r + halfWidth) * blurFilterWidth + c + halfWidth] = filterValue;
filterSum += filterValue;
}
}
float normalizationFactor = 1.f / filterSum;
for (int r = -halfWidth; r <= halfWidth; ++r)
{
for (int c = -halfWidth; c <= halfWidth; ++c)
{
(*filter)[(r + halfWidth) * blurFilterWidth + c + halfWidth] *= normalizationFactor;
}
}
}
void freeFilter(float *filter)
{
delete[] filter;
}
int main( int argc, char** argv )
{
using namespace cv;
using namespace std;
using namespace std::chrono;
if( argc != 2)
{
cout <<" Usage: blur_image imagefile" << endl;
return -1;
}
Mat image, blurredImage, referenceBlurredImage;
uchar4 *imageArray, *blurredImageArray;
prepareImagePointers(argv[1], image, &imageArray, blurredImage, &blurredImageArray, CV_8UC4);
int numRows = image.rows, numCols = image.cols;
float *filter, filterSigma;
int filterWidth;
prepareFilter(&filter, &filterWidth, &filterSigma);
cv::Size filterSize(filterWidth, filterWidth);
auto start = system_clock::now();
cv::GaussianBlur(image, referenceBlurredImage, filterSize, filterSigma, filterSigma, BORDER_REPLICATE);
auto duration = duration_cast<milliseconds>(system_clock::now() - start);
cout<<"OpenCV time (ms):" << duration.count() << endl;
start = system_clock::now();
BlurImageOpenMP(imageArray, blurredImageArray, numRows, numCols, filter, filterWidth);
duration = duration_cast<milliseconds>(system_clock::now() - start);
cout<<"OpenMP time (ms):" << duration.count() << endl;
cout<<"OpenMP similarity:" << getEuclidianSimilarity(referenceBlurredImage, blurredImage) << endl;
for (int i=0; i<4; ++i)
{
memset(blurredImageArray, 0, sizeof(uchar4)*numRows*numCols);
start = system_clock::now();
BlurImageCUDA(imageArray, blurredImageArray, numRows, numCols, filter, filterWidth);
duration = duration_cast<milliseconds>(system_clock::now() - start);
cout<<"CUDA time full (ms):" << duration.count() << endl;
cout<<"CUDA similarity:" << getEuclidianSimilarity(referenceBlurredImage, blurredImage) << endl;
}
freeFilter(filter);
return 0;
}
По пунктам:
double getEuclidianSimilarity(const cv::Mat& a, const cv::Mat& b)
{
double errorL2 = cv::norm(a, b, cv::NORM_L2);
double similarity = errorL2 / (double) (a.rows * a.cols);
return similarity;
}
По сути, она находит среднюю сумму квадратов разниц значений всех каналов всех пикселей двух изображений.
OpenMP реализация алгоритма:
#include <stdio.h>
#include <omp.h>
#include <vector_types.h>
#include <vector_functions.h>
void BlurImageOpenMP(const uchar4 * const imageArray,
uchar4 * const blurredImageArray,
const long numRows,
const long numCols,
const float * const filter,
const size_t filterWidth)
{
using namespace std;
const long halfWidth = filterWidth/2;
#pragma omp parallel for collapse(2)
for (long row = 0; row < numRows; ++row)
{
for (long col = 0; col < numCols; ++col)
{
float resR=0.0f, resG=0.0f, resB=0.0f;
for (long filterRow = -halfWidth; filterRow <= halfWidth; ++filterRow)
{
for (long filterCol = -halfWidth; filterCol <= halfWidth; ++filterCol)
{
//Find the global image position for this filter position
//clamp to boundary of the image
const long imageRow = min(max(row + filterRow, static_cast<long>(0)), numRows - 1);
const long imageCol = min(max(col + filterCol, static_cast<long>(0)), numCols - 1);
const uchar4 imagePixel = imageArray[imageRow*numCols+imageCol];
const float filterValue = filter[(filterRow+halfWidth)*filterWidth+filterCol+halfWidth];
resR += imagePixel.x*filterValue;
resG += imagePixel.y*filterValue;
resB += imagePixel.z*filterValue;
}
}
blurredImageArray[row*numCols+col] = make_uchar4(resR, resG, resB, 255);
}
}
}
Для всех 3 каналов каждого пиксела исходного изображения считаем описанную взвешенную сумму, результат записываем в соответствующую позицию выходного изображения.
CUDA вариант:
#include <iostream>
#include <cuda.h>
#include "CUDA_wrappers.hpp"
#include "common/CUDA_common.hpp"
__global__
void gaussian_blur(const uchar4* const d_image,
uchar4* const d_blurredImage,
const int numRows,
const int numCols,
const float * const d_filter,
const int filterWidth)
{
const int row = blockIdx.y*blockDim.y+threadIdx.y;
const int col = blockIdx.x*blockDim.x+threadIdx.x;
if (col >= numCols || row >= numRows)
return;
const int halfWidth = filterWidth/2;
extern __shared__ float shared_filter[];
if (threadIdx.y < filterWidth && threadIdx.x < filterWidth)
{
const int filterOff = threadIdx.y*filterWidth+threadIdx.x;
shared_filter[filterOff] = d_filter[filterOff];
}
__syncthreads();
float resR=0.0f, resG=0.0f, resB=0.0f;
for (int filterRow = -halfWidth; filterRow <= halfWidth; ++filterRow)
{
for (int filterCol = -halfWidth; filterCol <= halfWidth; ++filterCol)
{
//Find the global image position for this filter position
//clamp to boundary of the image
const int imageRow = min(max(row + filterRow, 0), numRows - 1);
const int imageCol = min(max(col + filterCol, 0), numCols - 1);
const uchar4 imagePixel = d_image[imageRow*numCols+imageCol];
const float filterValue = shared_filter[(filterRow+halfWidth)*filterWidth+filterCol+halfWidth];
resR += imagePixel.x * filterValue;
resG += imagePixel.y * filterValue;
resB += imagePixel.z * filterValue;
}
}
d_blurredImage[row*numCols+col] = make_uchar4(resR, resG, resB, 255);
}
void BlurImageCUDA(const uchar4 * const h_image,
uchar4 * const h_blurredImage,
const size_t numRows,
const size_t numCols,
const float * const h_filter,
const size_t filterWidth)
{
uchar4 *d_image, *d_blurredImage;
cudaSetDevice(0);
checkCudaErrors(cudaGetLastError());
const size_t numPixels = numRows * numCols;
const size_t imageSize = sizeof(uchar4) * numPixels;
//allocate memory on the device for both input and output
checkCudaErrors(cudaMalloc(&d_image, imageSize));
checkCudaErrors(cudaMalloc(&d_blurredImage, imageSize));
//copy input array to the GPU
checkCudaErrors(cudaMemcpy(d_image, h_image, imageSize, cudaMemcpyHostToDevice));
float *d_filter;
const size_t filterSize = sizeof(float) * filterWidth * filterWidth;
checkCudaErrors(cudaMalloc(&d_filter, filterSize));
checkCudaErrors(cudaMemcpy(d_filter, h_filter, filterSize, cudaMemcpyHostToDevice));
dim3 blockSize;
dim3 gridSize;
int threadNum;
cudaEvent_t start, stop;
cudaEventCreate(&start);
cudaEventCreate(&stop);
threadNum = 32;
blockSize = dim3(threadNum, threadNum, 1);
gridSize = dim3(numCols/threadNum+1, numRows/threadNum+1, 1);
cudaEventRecord(start);
gaussian_blur<<<gridSize, blockSize, filterSize>>>(d_image,
d_blurredImage,
numRows,
numCols,
d_filter,
filterWidth);
cudaEventRecord(stop);
cudaEventSynchronize(stop);
cudaDeviceSynchronize(); checkCudaErrors(cudaGetLastError());
float milliseconds = 0;
cudaEventElapsedTime(&milliseconds, start, stop);
std::cout << "CUDA time kernel (ms): " << milliseconds << std::endl;
checkCudaErrors(cudaMemcpy(h_blurredImage, d_blurredImage, sizeof(uchar4) * numPixels, cudaMemcpyDeviceToHost));
checkCudaErrors(cudaFree(d_filter));
checkCudaErrors(cudaFree(d_image));
checkCudaErrors(cudaFree(d_blurredImage));
}
extern __shared__ float shared_filter[];
if (threadIdx.y < filterWidth && threadIdx.x < filterWidth)
{
const int filterOff = threadIdx.y*filterWidth+threadIdx.x;
shared_filter[filterOff] = d_filter[filterOff];
}
__syncthreads();
Здесь __shared__ в объявлении массива весов фильтра говорит, что эти данные следует разместить в общей памяти блока; extern значит, что размер выделяемой памяти будет задан в месте вызова ядра; __syncthreads — барьер, гарантирующий что все веса фильтра будут перенесены в общую память блока прежде чем какой либо из потоков этого блока продолжит свое выполнение. В дальнейшем, все чтения весов фильтра осуществляются уже из общей памяти блока.
Компилируем, запускаем:
OpenCV time (ms):2
OpenMP time (ms):11
OpenMP similarity:0.00287131
CUDA time kernel (ms): 2.93245
CUDA time full (ms):32
CUDA similarity:0.00287131
CUDA time kernel (ms): 2.93402
CUDA time full (ms):4
CUDA similarity:0.00287131
CUDA time kernel (ms): 2.93267
CUDA time full (ms):4
CUDA similarity:0.00287131
CUDA time kernel (ms): 2.93312
CUDA time full (ms):4
CUDA similarity:0.00287131
Как видим, если не учитывать самый первый запуск CUDA-варианта, получили почти 3-х кратный выигрыш по времени в сравнении с OpenMP вариантом, хотя и не догнали OpenCV вариант — который, кстати, использует OpenCL [5].
На сегодня все, в следующей части рассмотрим некоторые фундаментальные алгоритмы GPU.
Весь исходный код доступен на bitbucket [6].
Автор: VladGorbatiuk
Источник [7]
Сайт-источник PVSM.RU: https://www.pvsm.ru
Путь до страницы источника: https://www.pvsm.ru/udacity/76863
Ссылки в тексте:
[1] Часть 1: Введение.: http://habrahabr.ru/company/epam_systems/blog/245503/
[2] Fermi: http://en.wikipedia.org/wiki/Fermi_(microarchitecture)
[3] тут: http://docs.nvidia.com/cuda/cuda-c-programming-guide/#atomic-functions
[4] по Гауссу: http://en.wikipedia.org/wiki/Gaussian_blur
[5] OpenCL: http://en.wikipedia.org/wiki/OpenCL
[6] bitbucket: https://bitbucket.org/VladyslavGorbatiuk/cuda-parallel-programming-code
[7] Источник: http://habrahabr.ru/post/245523/
Нажмите здесь для печати.