М. В. Ломоносова 2015 Антонюк Валерий Алексеевич Программирование на видеокартах (gpgpu). Спецкурс кафедры мми. М. Физический факультет мгу им. М. В. Ломоносова, 2015. 48 с. Учебное пособие
 Скачать 1.54 Mb.
|
|
Умножение матриц Примером весьма часто используемого вычислительного алгоритма (но не очень хорошо поддающегося распараллеливанию) является алгоритм умножения матриц. Стандартный вид алгоритма для случая выполнения на одном процессоре таков: for(i = 0; i < M; i++) for(j = 0; j < N; j++) for(k = 0; k < K; k++) C[i][j] += A[i][k] * B[k][j]; Таким образом, результирующая матрица формируется поэлементно (здесь порядок расчёта элементов — вследствие коммутативности используемых операций — значения не имеет, по крайней мере — теоретически), каждый элемент вычисляется на основе одной строки левой матрицы и одного столбца правой матрицы. Сами матрицы — для простоты — в этом иллюстративном примере адресуются как двумерные массивы, хотя на практике элементы матриц хранятся (чаще всего — по строкам) в одномерных массивах, указатели на которые вместе с размерами матриц образуют структуры вроде такой: typedef struct { int width; int height; float* elements; } Matrix; Ниже приведена функция-ядро для реализации подобного «наивного» способа перемножения матриц ( C = A*B ); предполагается, что матрицы расположены в глобальной памяти. Каждая нить считывает строку row из первой матрицы и столбец col из второй и вычисляет «свой» элемент матрицы ( row,col ): __global__ void MatMulKernel(Matrix A, Matrix B, Matrix C){ // каждая нить вычисляет один элемент матрицы C float Cvalue = 0; int row = blockIdx.y * blockDim.y + threadIdx.y; int col = blockIdx.x * blockDim.x + threadIdx.x; for (int e = 0; e < A.width; ++e) Cvalue += A.elements[row*A.width + e] * B.elements[e*B.width + col]; C.elements[row*C.width + col] = Cvalue; } Подобная реализация выглядит привычно и несложно, однако, имеет один существенный недостаток: для формирования каждого элемента результирующей матрицы из глобальной памяти многократно читаются (и перечитываются) элементы отдельных строк и столбцов. Для того, чтобы уменьшить число загрузок элементов исходных матриц, имеет смысл это делать порциями (скажем, квадратными блоками), которые могут быть временно — для формирования произведения этих блоков — помещены из глобальной памяти в более быструю разделяемую память. А поскольку теперь надо будет работать не только с исходными матрицами, но и блоками этих матриц, структуру описания отдельной матрицы следует дополнить ещё одним параметром, содержащим ширину исходной матрицы (в приводимом далее тексте это поле структуры именуется stride ). 3 14 Иллюстрация алгоритма перемножения матриц из руководства NVIDIA. При использовании разделяемой ( __shared__ ) памяти каждый блок нитей «занимается» вычислением квадратной подматрицы Csub (размера BLOCKSIZE на BLOCKSIZE ), а каждая нить в этом блоке вычисляет один элемент подматрицы Csub . Подматрица Csub есть произведение двух прямоугольных матриц: подматрицы матрицы A (с размерами A.width на BLOCKSIZE ) и подматрицы матрицы B (с размерами BLOCKSIZE на B.height ). Из-за ограниченного размера разделяемой памяти эти прямоугольные матрицы приходится разбивать на необходимое количество квадратных подматриц размера BLOCKSIZE и подматрица Csub вычисляется как сумма произведений этих квадратных матриц (это возможно как следствие правила блочного перемножения матриц). Каждое из таких произведений вычисляется после загрузки двух соответствующих квадратных матриц из глобальной памяти в разделяемую, причём сначала каждая нить загружает один элемент каждой матрицы, а затем вычисляет один элемент их произведения. Отдельные нити аккумулируют результат каждого из произведений в регистре, а по завершении записывают этот результат в глобальную память. В реализации алгоритма используются три вспомогательные функции: GetElement() , SetElement() и GetSubMatrix() , предназначенные для работы в GPU (они помечены спецификатором __device__ и по этой причине могут быть вызваны только из ядра). Надо сказать, что хотя — для удобства программирования — эти фрагменты кода и оформлены как отдельные функции, на самом деле они просто «вставляются» в код, поскольку вызовы в системе команд GPU обычно отсутствуют. Первая функция обеспечивает получение элемента из матрицы или подматрицы: __device__ float GetElement(const Matrix A, int row, int col){ return A.elements[row*A.stride + col]; } Вторая функция заносит заданное значение в элемент матрицы или подматрицы: __device__ void SetElement(Matrix A, int row, int col, float value){ A.elements[row*A.stride + col] = value; } 15 Третья функция реализует извлечение подматрицы-блока из исходной матрицы: __device__ Matrix GetSubMatrix(Matrix A, int row, int col){ Matrix ASub; ASub.width = BLOCK_SIZE; ASub.height = BLOCK_SIZE; ASub.stride = A.stride; ASub.elements = &A.elements[A.stride*BLOCK_SIZE*row + BLOCK_SIZE*col]; return ASub; } А вот, собственно, и сам алгоритм такого блочного перемножения матриц (размеры которых предполагаются кратными размерам блока): __global__ void MatMulKernel(Matrix A, Matrix B, Matrix C){ // Координаты блока в матрице ( номер блока в строке и в столбце ) int blockRow = blockIdx.y; int blockCol = blockIdx.x; // Каждый блок нитей вычисляет одну подматрицу Csub, при этом // каждая нить создаёт свой собственный дескриптор матрицы Csub Matrix Csub = GetSubMatrix(C, blockRow, blockCol); // Каждая нить вычисляет один элемент подматрицы Csub float Cvalue = 0; // thread row and col WITHIN CSUB int row = threadIdx.y; int col = threadIdx.x; // Цикл по всем подматрицам строки блоков A и столбца блоков B; // этот цикл необходим для вычисления всей подматрицы Csub. // Блочное умножение пары подматриц и аккумуляция результата for (int m = 0; m < (A.width / BLOCK_SIZE); ++m){ // Формирование дескрипторов подматриц Asub и Bsub Matrix Asub = GetSubMatrix(A, blockRow, m); Matrix Bsub = GetSubMatrix(B, m, blockCol); // Копирование элементов ASub и Bsub в разделяемую память // Каждая нить загружает один элемент ASub и один – Bsub // Обратите внимание : каждая нить объявляет матрицы As и Bs, // хотя блок нитей содержит только одну матрицу As и одну Bs __shared__ float As[BLOCK_SIZE][BLOCK_SIZE]; __shared__ float Bs[BLOCK_SIZE][BLOCK_SIZE]; As[row][col] = GetElement(Asub, row, col); Bs[row][col] = GetElement(Bsub, row, col); // Нити синхронизируются – чтобы убедиться , что всё прочитано __syncthreads(); // Скалярное произведение одной строки Asub и одного столбца Bsub // формируют ( частично ) один элемент результирующей подматрицы for(int e=0; e // Надо убедиться , что все Cvalues просуммированы до начала // считывания последующих блоков в разделяемые матрицы As, Bs __syncthreads(); } // Запись Csub в глобальную память : каждая нить записывает « свой » элемент SetElement(Csub, row , col, Cvalue); } 16 Алгоритм параллельной сортировки Bitonic Sort Слово bitonic относится к последовательности сортируемых числовых величин и обозначает такую последовательность, которая сначала монотонно возрастает, а затем монотонно убывает, или наоборот (либо циклически приводимую к ней). В дальнейшем предполагаем, что количество величин является степенью двойки. Если обменивать между собой соответствующие величины из первой и второй половин (в случае их неправильного положения), мы придём к ситуации, где битоническими будут обе половины. Это так называемый шаг битонического расщепления. Проделаем аналогичные действия для каждой из половин, затем для половин этих половин и т.д. Обратите внимание, что на самом первом шаге в первой половине удалось сосредоточить все величины, которые уже больше не придётся обменивать со второй половиной! Таким образом, битоническая сортировка использует свойство битонического расщепления: величины отсортированы для половин, а каждая половина является битонической. Повторно применяя битоническое расщепление мы можем превратить битоническую последовательность в монотонную (т.е., отсортированную), когда в каждой соседней паре величины располагаются в "правильном" порядке, а сами пары уже расположены "правильно". Для сортировки произвольной последовательности её сначала надо превратить в битоническую, а затем применять битоническое расщепление. Осуществляется это путём поочерёдного сравнения и перестановки сначала соседних пар (с чередованием порядка следования), затем четвёрок (сначала — удалённых пар, потом — соседних), затем аналогичным образом восьмёрок и так далее. Пары, используемые для сравнения и перестановки на каждом шаге, показаны на рисунке: Иллюстрация действий в Bitonic Sort ( http://en.wikipedia.org/wiki/Bitonic_sorter ) Оказывается, все эти действия укладываются в однотипную последовательность шагов, описываемую двумя вложенными циклами (см. код далее). Параметры для Bitonic Sort с помощью CUDA: /* Каждая нить " работает " с одной величиной */ #define THREADS 512 // 2^9 #define BLOCKS 32768 // 2^15 #define NUM_VALS THREADS*BLOCKS 4 17 Функция-ядро: __global__ void bitonic_sort_step(float *dev_values, int j, int k) { unsigned int i, ixj; /* Сортируемые " партнёры ": i and ixj */ i = threadIdx.x + blockDim.x * blockIdx.x; ixj = i^j; if ((ixj)>i) { if ((i&k)==0) { /* Сортировка по возрастанию */ if (dev_values[i]>dev_values[ixj]) { /* Обмен i и ixj */ float temp = dev_values[i]; dev_values[i] = dev_values[ixj]; dev_values[ixj] = temp; } } if ((i&k)!=0) { /* Сортировка по убыванию */ if (dev_values[i] Обмен i и ixj */ float temp = dev_values[i]; dev_values[i] = dev_values[ixj]; dev_values[ixj] = temp; } } } } А вот сама процедура сортировки, включающая формирование битонической последовательности из произвольной заданной: /** * Bitonic Sort - сортировка происходит " на месте "! */ void bitonic_sort(float *values) { float *dev_values; size_t size = NUM_VALS * sizeof(float); cudaMalloc((void**) &dev_values, size); cudaMemcpy(dev_values, values, size, cudaMemcpyHostToDevice); dim3 blocks(BLOCKS,1); /* Number of blocks */ dim3 threads(THREADS,1); /* Number of threads */ int j, k; /* k - " длина " первого шага сравнения - обмена ; он удваивается */ for (k = 2; k <= NUM_VALS; k <<= 1) { 18 /* Перестановки должны продолжаться пока шаг не уменьшился */ for (j=k>>1; j>0; j=j>>1) { bitonic_sort_step<< } } cudaMemcpy(values, dev_values, size, cudaMemcpyDeviceToHost); cudaFree(dev_values); } Выделяется память для сортируемых величин на видеокарте (GPU) и они копируются туда. Затем несколько раз вызывается функция-ядро для различных параметров j, k . И всё оказывается отсортированным, причём прямо "на месте"! Результат копируется с видеокарты на компьютер, выделенная ранее память на карте — освобождается. "Косвенное" использование возможностей GPU Не всегда у программирующего есть возможность (и желание) изучать "прямое" программирование графической карты. Это может быть связано с изменчивостью архитектуры и частой сменой поколений графических карт, завязанностью такого выбора на одного конкретного производителя карт, сложностью самой архитектуры. К счастью, довольно часто можно не углубляться в изучение "железок", а воспользоваться тем, что наиболее употребительные алгоритмы уже реализованы в рамках либо поставляемых с CUDA, либо каких-то сторонних библиотек. В состав CUDA входят такие специализированные библиотеки, как Thrust, cuBLAS, cuFFT, cuRAND и т.п. Библиотека Thrust — это средства работы с контейнером vector на графической карте (GPU) в стиле библиотеки STL. Здесь можно не отвлекаться на выделение/освобождение памяти на компьютере (CPU) или на графической карте (GPU), а просто объявлять необходимые вектора, причём освобождаться они будут автоматически. Выделение памяти для вектора в памяти компьютера: thrust::host_vector Генерация случайных значений (в стиле библиотеки STL): thrust::generate(h_vec.begin(), h_vec.end(), rand); Создание места для вектора на графической карте и его копирование: thrust::device_vector Сортировка (на видеокарте; прямо в месте расположения вектора): thrust::sort(d_vec.begin(), d_vec.end()); Обратное копирование результата с видеокарты в память компьютера: thrust::copy(d_vec.begin(), d_vec.end(), h_vec.begin()); Видно, что здесь используется явное указание пространства имён, поскольку названия совпадают с аналогичными названиями из библиотеки STL. 19 Ещё один пример: умножение вектора на число и суммирование с другим вектором. В стиле CUDA C: __global__ void saxpy_kernel(int n, float a, float *x, float *y) { const int i = blockDim.x*blockIdx.x+threadIdx.x; if (i < n) y[i] = a*x[i]+y[i]; } void saxpy(int n, float a, float *x, float *y) { // параметры запуска ядра int block_size = 256; int grid_size = (n + block_size-1)/block_size; // запуск ядра saxpy_kernel saxpy_kernel<< } В стиле Thrust: struct saxpy_functor { const float a; saxpy_functor(float _a) : a(_a) {} __host__ __device__ float operator()(float x, float y) { return a*x+y; } }; void saxpy(float a, device_vector { // определение функтора saxpy_functor func(a); // вызов преобразования transform(x.begin(), x.end(), y.begin(), y.begin(), func); } Напомним, что функтор — это переменная, к которой можно приписать круглые скобки со списком параметров в них, что приведёт к вызову оператора, обозначаемого круглыми скобками ( operator() ). Может быть также именем класса (или структуры). Если функтор — шаблонный, то имя шаблона с типом тоже будет функтором. Обратите внимание, что оператор "круглые скобки" помечен спецификаторами и __host__ , и __device__ , что означает, что он может применяться как в host-векторах, так и в device- векторах. Вывод результирующего вектора V (на компьютере) легко делается с помощью STL: std::copy(V.begin(), V.end(), std::ostream_iterator Библиотеки, использующие CUDA GPU Попытаемся вкратце охарактеризовать вычислительные библиотеки, использующие CUDA-карты. Некоторые из таких библиотек в настоящее время уже являются частью CUDA, поэтому имеет смысл сначала поговорить о них. Разберёмся сначала с Thrust, cuBLAS и cuSPARSE. Thrust: состав и возможности Библиотека Thrust является библиотекой шаблонов C++ для CUDA, аналогичной STL, но пока более простой. Она содержит только один вид контейнеров — векторы, однако два типа их: thrust::host_vector и thrust::device_vector . Первый (как следует из названия) предназначен для работы в памяти компьютера, второй — на графической карте. Подобно вектору в STL — это общие контейнеры, хранящие любые типы данных и избавляющие программиста от забот с выделением/освобождением памяти, а также упрощающие обмен данными между CPU и GPU. С помощью Thrust описываются собственно вычисления, а не то, как будут храниться данные или производиться эти вычисления. В библиотеке обеспечивается абстрактный интерфейс к фундаментальным параллельным алгоритмам, таким как сортировка ( thrust::sort() ), редуцирование ( thrust::reduce() ) и др. Широко используется принятое в STL указание диапазона как пары итераторов. Сами итераторы могут использоваться и как указатели, в том числе на память в GPU (например, для передачи в функцию-ядро): // выделяем память на графическом устройстве thrust::device_vector // получаем указатель на этот вектор в памяти GPU int * ptr = thrust::raw_pointer_cast(&d_vec[0]); // используем указатель в функции - ядре my_kernel<< // Операция разыменования в памяти CPU не имеет смысла ! Видно, что библиотека использует собственное пространство имён ( thrust ), что позволяет избежать коллизий с аналогичными именами |