Программирование на видеокартах GPGPU (1184391), страница 3
Текст из файла (страница 3)
В простейших случаях, когда достаточно одномернойорганизации данных, можно обойтись и двумя целочисленными значениями.9Более сложные случаи запуска могут иметь также и дополнительные параметры; вот как выглядитполный синтаксис запуска ядер на исполнение:kernelName<<<Dg,Db,Ns,S>>>(args);Здесь kernelName — имя __global__-функции, Dg — переменная типа dim3 (параметрысетки блоков), Db — переменная типа dim3 (параметры блока нитей), Ns — дополнительныйобъём разделяемой памяти в байтах для каждого блока (необязательный параметр, он поумолчанию равен нулю), S — задаёт поток (CUDAstream), в котором должен произойти вызов(поток 0 по умолчанию), args — параметры функции-ядра (их может быть несколько, а общийразмер этих параметров пока ограничен 256 байтами).Дополнительные математические функцииПоскольку ядра исполняются на GPU, а не на CPU, им необходимы альтернативные реализациикак минимум всех тех математических функций, которые были доступны в программах ранее врамках стандартной математической библиотеки.
Кроме double-версий в CUDA также имеютсяfloat-аналоги стандартных функций (например, для sin — sinf). Есть ещё набор функцийпониженной точности, но более «быстрых» (для sinf, например, __sinf).Для функций, реализующих арифметические операции и преобразование типов, можно такжезадавать с помощью суффиксов их имени способ округления результата: rn (к ближайшему),rz (к нулю), ru (вверх), rd (вниз).Разделяемая (__shared__) памятьОсновная часть динамической памяти GPU доступна как глобальная __device__-память — всемнитям и блокам, а также CPU.
Несмотря на то, что пересылка из памяти CPU в память GPUосуществляется довольно быстро, пересылки в рамках динамической памяти ещё быстрее.Разделяемая (__shared__) память располагается в GPU и потому обладает гораздо большимбыстродействием, чем глобальная, но доступна для чтения и записи нитям из одного блока.Первый способ выделения разделяемой памяти в ядре — явное указание размера массива соспецификатором __shared__.Второй способ — при запуске ядра задать дополнительный объём разделяемой памяти в байтах,который необходимо выделить каждому блоку при запуске ядра (третий параметр конфигурациифункции-ядра, см. выше). Для доступа в ядрах к такой памяти используется описание массива вфункции ядра без явно заданного размера, например:__shared__ float buf[];При этом предполагается, что сам вызов ядра (с названием kernel) выглядит так:kernel<<< ... , ...
, k*sizeof(float) >>>( ... );Тогда каждый блок ядер будет иметь доступ к массиву из k вещественных значений вразделяемой памяти (см. реализацию алгоритма умножения матриц далее).10Взаимодействие нитей в блокахНити в рамках блока могут взаимодействовать между собой с помощью общей памяти, а потомупоявляется необходимость синхронизировать выполнение кода — чтобы чтение результатов изпамяти не происходило ранее их формирования. Для этого программист определяет в коде точкисинхронизации, используя вызов функции __syncthreads(). Она действует как барьер, передкоторым все нити блока должны остановиться, ожидая завершения работы остальных нитейблока.Адресация памяти в нитяхИспользуемая в NVIDIA GPU архитектура относится к типу SIMT (Single Instruction MultipleThreads, «одна инструкция, много нитей»): один и тот же код исполняется большимколичеством нитей. Для эффективного распараллеливания таких вычислений, где данныемогут обрабатываться одинаково и независимо, нужно, чтобы каждая нить обрабатывала"свою" часть данных, поэтому при программной реализации подобных вычислений возникаетнеобходимость "правильного" распределения данных по нитям.
Добиться этого можно такойпроцедурой сопоставления нитей и участков памяти для них, где разным нитям соответствуютразные участки памяти. При этом часто удобно, чтобы нитям одного блока сопоставлялся одинкомпактный кусок памяти (без пропусков).Как оказывается, каждая копия функции ядра после запуска на исполнение имеет доступ кспециальным переменным: к размерностям блоков и сетки (они называются blockDim иgridDim соответственно, обе имеют тип dim3), а также к положениям конкретной копии ядрав блоке и этого блока во всей сетке при исполнении кода (threadIdx, blockIdx, обепеременные имеют тип int3).
Используя эти переменные, каждая копия ядра можетсформировать уникальные индексы для доступа к данным.Можно, например, последовательно перебирать (т.е. нумеровать) нити в блоке по каждомувозможному измерению по очереди, при этом количество использованных индексов по одномуизмерению будет равно размеру блока по этому измерению, а общее количество использованныхномеров-индексов будет равно произведению размеров блока по всем возможным измерениям.Аналогично можно поступить и с блоками в сетке.Порядок перебора измерений при этом часто не важен, хотя традиционно всё начинается скомпоненты .x, затем — если необходимо — используется компонента .y, потом — еслинадо — используется компонента .z. Объяснение этому простое: разные измерения вконфигурации запуска ядер, как правило, неравноправны по предельным возможностям,поэтому «традиционный» выбор порядка следования измерений может упрощатьпоследующее масштабирование программы.Подобная схема адресации часто реализуется в программах.
Рассмотрим несколько примеровформирования индексов для доступа к данным в зависимости от объявленной конфигурацииисполнения ядер.Например, в программе bitreverse.cu, осуществляющей "перевёртывание" битов в каждомиз байтов последовательности (организованной в четвёрки целых "слов"), конфигурациясодержит всего один блок с одномерной организацией нитей, поэтому для доступа к "словам"достаточно одной компоненты индекса нити (threadIdx.x). В программе pi.cu,вычисляющей приближение к числу π суммированием ряда, конфигурация нитей тоже11одномерна, равно как и сетка блоков, поэтому для доступа к местам подсчёта отдельныхслагаемых ряда используется уникальный индекс, составленный из индекса нити поединственному измерению в блоке (threadIdx.x) и индекса блока по единственномуизмерению в сетке (blockIdx.x):int idx = blockIdx.x*blockDim.x+threadIdx.x;Выбор blockDim.x в качестве сомножителя в этом выражении обеспечивает "неразрывность"изменения уникального индекса (threadIdx.x изменяется от 0 до blockDim.x-1).В программе MatrixAdd.cu, предназначенной для суммирования квадратных матриц NxN,конфигурация запуска ядер двумерна и по нитям в блоках, и по блокам в сетке — сообразносмыслу задачи, — а сама матрица разбита на квадратные подблоки такого размера (16x16), чтобыкаждый подблок обрабатывался одним блоком нитей.
Поэтому формирование уникальныхиндексов здесь организовано по-другому: создаются индексы сначала по каждому измерению (наоснове нужных компонент индекса нити в блоке и индекса блока в сетке), а затем оникомбинируются в один общий индекс для адресации элементов матрицы, расположенной влинейной памяти.int i = blockIdx.x * blockDim.x + threadIdx.x;int j = blockIdx.y * blockDim.y + threadIdx.y;int index = i + j*N;В программе compute_pi.cu (ещё один вариант вычисления частичной суммы ряда для π)используется трёхмерная конфигурация нитей в единственном блоке, поэтому формула дляуникального индекса использует только компоненты положения нити в блоке, причём переборизмерений здесь производится немного непривычным (хотя и вполне допустимым) образом, т.е.,сначала — компонента .z, потом — .y, а затем только — .x:myid = threadIdx.z + threadDim*threadIdx.y+threadDim*threadIdx.y*threadIdx.x;Теперь уже понятно, что общая формула адресации нити в блоке (получение порядкового индексанити tid) может быть такова:tid = threadIdx.x +threadIdx.y * blockDim.x +threadIdx.z * blockDim.x * blockDim.yОна годится для любой возможной размерности конфигурации нитей в блоке, потому что, есликонфигурация двумерна (blockDim.z = 1), то единственным возможным значением дляthreadIdx.z будет 0; если же она одномерна (т.е., ещё и blockDim.y = 1), то иthreadIdx.y будет всегда равно 0, при этом реальный вид формулы для двумерной илиодномерной организации нитей будет проще (и совпадёт с формулами в рассмотренныхпримерах программ).Аналогично, формула адресации блока в сетке (получение порядкового индекса блока bid)может быть подобной:bid = blockIdx.x +blockIdx.y * gridDim.x +blockIdx.z * gridDim.x * gridDim.y12Она тоже годится для любой возможной размерности конфигурации блоков в сетке (еслиgridDim.z = 1, то blockIdx.z = 0; если же ещё и gridDim.y = 1, то иblockIdx.y = 0).Вводя вспомогательную величину nthreads для общего числа нитей в каждом из блоков:nthreads = blockDim.x * blockDim.y * blockDim.zмы можем получить линейный индекс id (поскольку память адресуется линейно) для доступа кданным из произвольной нити произвольного блока:id = tid + nthreads * bidНо иногда (как показывает пример с матрицей) удобно отступить от вышеизложенной схемы,организуя как "сквозную" адресацию нитей по каждому измерению, так и "блочную" — чтобыиндексы нитей по каждому из двух измерений соответствовали индексам элементов матрицыв отдельном матричном блоке, индексы блоков нитей соответствовали индексам блоковматрицы, а формируемые по каждому измерению комбинированные индексы — индексамэлементов матрицы в целом.Суммирование компонент вектора — пример ядраРассмотрим в качестве примера один из вариантов реализации суммирования компонент вектора(обобщённо реализуемая операция называется в англоязычной литературе reduction и можетиспользовать не только сложение, но и другие бинарные коммутативные операции, например,умножение, логическое сложение и т.п., приводя к формированию произведения и т.д.):__global__ void plus_reduce(int *g_idata, int N, int *g_odata) {__shared__ int sdata[BLOCKSIZE];// Каждая нить загружает одно значение из глобальной памятиunsigned int tid = threadIdx.x;unsigned int i = blockIdx.x*blockDim.x + threadIdx.x;sdata[tid] = (i < N)? g_idata[i] : 0;__syncthreads();// Суммирование осуществляется в разделяемой памятиfor (unsigned int s=blockDim.x/2; s>0; s>>=1) {if (tid < s)sdata[tid] += sdata[tid + s];__syncthreads();}// Запись результата для данного блока в глобальную памятьif (tid == 0) g_odata[blockIdx.x] = sdata[0];}Здесь предполагается, что окончательное получение суммы осуществляется CPU.