Программирование на видеокартах GPGPU (1184391), страница 4
Текст из файла (страница 4)
Второйвариант: добавление сразу к общей сумме, находящейся в первом элементе.if (tid == 0) atomicAdd(g_odata, sdata[0]);Т.е., может использоваться одна из так называемых атомарных операций, гарантирующаякорректность выполнения операции в том случае, когда её пытаются выполнить многие нитиодновременно.133Умножение матрицПримером весьма часто используемого вычислительного алгоритма (но не очень хорошоподдающегося распараллеливанию) является алгоритм умножения матриц.Стандартный вид алгоритма для случая выполнения на одном процессоре таков: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){// каждая нить вычисляет один элемент матрицы Cfloat 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).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, при этом// каждая нить создаёт свой собственный дескриптор матрицы CsubMatrix Csub = GetSubMatrix(C, blockRow, blockCol);// Каждая нить вычисляет один элемент подматрицы Csubfloat Cvalue = 0;// thread row and col WITHIN CSUBint row = threadIdx.y;int col = threadIdx.x;// Цикл по всем подматрицам строки блоков A и столбца блоков B;// этот цикл необходим для вычисления всей подматрицы Csub.// Блочное умножение пары подматриц и аккумуляция результатаfor (int m = 0; m < (A.width / BLOCK_SIZE); ++m){// Формирование дескрипторов подматриц Asub и BsubMatrix 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();// Скалярное// формируютfor(int e=0;Cvalue +=произведение одной строки Asub и одного столбца Bsub(частично) один элемент результирующей подматрицыe<BLOCK_SIZE; ++e)As[row][e] * Bs[e][col];// Надо убедиться, что все Cvalues просуммированы до начала// считывания последующих блоков в разделяемые матрицы As, Bs__syncthreads();}// Запись Csub в глобальную память: каждая нить записывает «свой» элементSetElement(Csub, row , col, Cvalue);}164Алгоритм параллельной сортировки Bitonic SortСлово bitonic относится к последовательности сортируемых числовых величин и обозначает такуюпоследовательность, которая сначала монотонно возрастает, а затем монотонно убывает, илинаоборот (либо циклически приводимую к ней).
В дальнейшем предполагаем, что количествовеличин является степенью двойки.Если обменивать между собой соответствующие величины из первой и второй половин (в случаеих неправильного положения), мы придём к ситуации, где битоническими будут обе половины.Это так называемый шаг битонического расщепления. Проделаем аналогичные действия длякаждой из половин, затем для половин этих половин и т.д.
Обратите внимание, что на самомпервом шаге в первой половине удалось сосредоточить все величины, которые уже больше непридётся обменивать со второй половиной!Таким образом, битоническая сортировка использует свойство битонического расщепления:величины отсортированы для половин, а каждая половина является битонической. Повторноприменяя битоническое расщепление мы можем превратить битоническую последовательность вмонотонную (т.е., отсортированную), когда в каждой соседней паре величины располагаются в"правильном" порядке, а сами пары уже расположены "правильно".Для сортировки произвольной последовательности её сначала надо превратить в битоническую, азатем применять битоническое расщепление.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
