CUDA Reduction: Warp Unrolling (Школа)
В настоящее время я работаю над проектом, в котором развертываю последний перекос сокращения. Я закончил код выше; Тем не менее, некоторые изменения были сделаны путем угадывания, и я хотел бы объяснить, почему. Код, который я написал, является только функцией kernel4
// in is input array, out is where to store result, n is number of elements from in
// T is a float (32bit)
__global__ void kernel4(T *in, T *out, unsigned int n)
который является алгоритмом сокращения, остальная часть кода уже была предоставлена.
Код:
#include <stdlib.h>
#include <stdio.h>
#include "timer.h"
#include "cuda_utils.h"
typedef float T;
#define N_ (8 * 1024 * 1024)
#define MAX_THREADS 256
#define MAX_BLOCKS 64
#define MIN(x,y) ((x < y) ? x : y)
#define tid threadIdx.x
#define bid blockIdx.x
#define bdim blockDim.x
#define warp_size 32
unsigned int nextPow2( unsigned int x ) {
--x;
x |= x >> 1;
x |= x >> 2;
x |= x >> 4;
x |= x >> 8;
x |= x >> 16;
return ++x;
}
void getNumBlocksAndThreads(int whichKernel, int n, int maxBlocks, int maxThreads, int &blocks, int &threads)
{
if (whichKernel < 3) {
threads = (n < maxThreads) ? nextPow2(n) : maxThreads;
blocks = (n + threads - 1) / threads;
} else {
threads = (n < maxThreads*2) ? nextPow2((n + 1)/ 2) : maxThreads;
blocks = (n + (threads * 2 - 1)) / (threads * 2);
}
if (whichKernel == 5)
blocks = MIN(maxBlocks, blocks);
}
T reduce_cpu(T *data, int n) {
T sum = data[0];
T c = (T) 0.0;
for (int i = 1; i < n; i++)
{
T y = data[i] - c;
T t = sum + y;
c = (t - sum) - y;
sum = t;
}
return sum;
}
__global__ void
kernel4(T *in, T *out, unsigned int n)
{
__shared__ volatile T d[MAX_THREADS];
unsigned int i = bid * bdim + tid;
n >>= 1;
d[tid] = (i < n) ? in[i] + in[i+n] : 0;
__syncthreads ();
for(unsigned int s = bdim >> 1; s > warp_size; s >>= 1) {
if(tid < s)
d[tid] += d[tid + s];
__syncthreads ();
}
if (tid < warp_size) {
if (n > 64) d[tid] += d[tid + 32];
if (n > 32) d[tid] += d[tid + 16];
d[tid] += d[tid + 8];
d[tid] += d[tid + 4];
d[tid] += d[tid + 2];
d[tid] += d[tid + 1];
}
if(tid == 0)
out[bid] = d[0];
}
int main(int argc, char** argv)
{
T *h_idata, h_odata, h_cpu;
T *d_idata, *d_odata;
struct stopwatch_t* timer = NULL;
long double t_kernel_4, t_cpu;
int whichKernel = 4, threads, blocks, N, i;
if(argc > 1) {
N = atoi (argv[1]);
printf("N: %d\n", N);
} else {
N = N_;
printf("N: %d\n", N);
}
getNumBlocksAndThreads (whichKernel, N, MAX_BLOCKS, MAX_THREADS, blocks, threads);
stopwatch_init ();
timer = stopwatch_create ();
h_idata = (T*) malloc (N * sizeof (T));
CUDA_CHECK_ERROR (cudaMalloc (&d_idata, N * sizeof (T)));
CUDA_CHECK_ERROR (cudaMalloc (&d_odata, blocks * sizeof (T)));
srand48(time(NULL));
for(i = 0; i < N; i++)
h_idata[i] = drand48() / 100000;
CUDA_CHECK_ERROR (cudaMemcpy (d_idata, h_idata, N * sizeof (T), cudaMemcpyHostToDevice));
dim3 gb(blocks, 1, 1);
dim3 tb(threads, 1, 1);
kernel4 <<<gb, tb>>> (d_idata, d_odata, N);
cudaThreadSynchronize ();
stopwatch_start (timer);
kernel4 <<<gb, tb>>> (d_idata, d_odata, N);
int s = blocks;
while(s > 1) {
threads = 0;
blocks = 0;
getNumBlocksAndThreads (whichKernel, s, MAX_BLOCKS, MAX_THREADS, blocks, threads);
dim3 gb(blocks, 1, 1);
dim3 tb(threads, 1, 1);
kernel4 <<<gb, tb>>> (d_odata, d_odata, s);
s = (s + threads * 2 - 1) / (threads * 2);
}
cudaThreadSynchronize ();
t_kernel_4 = stopwatch_stop (timer);
fprintf (stdout, "Time to execute unrolled GPU reduction kernel: %Lg secs\n", t_kernel_4);
double bw = (N * sizeof(T)) / (t_kernel_4 * 1e9); // total bits / time
fprintf (stdout, "Effective bandwidth: %.2lf GB/s\n", bw);
CUDA_CHECK_ERROR (cudaMemcpy (&h_odata, d_odata, sizeof (T), cudaMemcpyDeviceToHost));
stopwatch_start (timer);
h_cpu = reduce_cpu (h_idata, N);
t_cpu = stopwatch_stop (timer);
fprintf (stdout, "Time to execute naive CPU reduction: %Lg secs\n", t_cpu);
if(abs (h_odata - h_cpu) > 1e-5)
fprintf(stderr, "FAILURE: GPU: %f CPU: %f\n", h_odata, h_cpu);
else
printf("SUCCESS: GPU: %f CPU: %f\n", h_odata, h_cpu);
return 0;
}
Мой первый вопрос: когда декларируешь
__shared__ volatile T d[MAX_THREADS];
Я хотел бы проверить мое понимание изменчивости. Volatile предотвращает неправильную оптимизацию компиляторами моего кода и обещает, что загрузка / сохранение завершаются через кеш, а не только через регистры (пожалуйста, исправьте меня, если ошиблись). Для сокращения, если суммы частичного сокращения все еще сохраняются в регистрах, почему это проблема?
Мой второй вопрос: когда делаю фактическое уменьшение деформации
if (tid < warp_size) { // Final log2(32) = 5 strides
if (n > 64) d[tid] += d[tid + 32];
if (n > 32) d[tid] += d[tid + 16];
d[tid] += d[tid + 8];
d[tid] += d[tid + 4];
d[tid] += d[tid + 2];
d[tid] += d[tid + 1];
}
Сумма редукции даст неверные результаты без условий (n > 64) и (n > 32). Результаты, которые я получаю:
FAILURE: GPU: 41.966557 CPU: 41.946209
При 5 испытаниях уменьшение графического процессора последовательно дает ошибку 0,0204. Я опасаюсь, что это ошибка операции с плавающей запятой.
Честно говоря, помощник моего учителя предложил это изменение, чтобы добавить условия (n > 64) и (n > 32), но не объяснил, почему это исправит код.
Поскольку n в моих исследованиях превышает 64, почему это обусловливает изменение результатов. У меня проблемы с поиском решения проблемы, потому что я не могу использовать функции печати, как в ЦП.
1 ответ
Давайте начнем с нескольких предварительных комментариев, прежде чем мы ответим на два ваших вопроса:
- Я рекомендую вам прочесть учебник NVIDIA по каноническому сокращению
- Сокращения, написанные таким образом, делают несколько предположений, одним из которых является то, что размер блока равен степени 2 (для "правильности").
- Ваш код использует варп-синхронное программирование на финальной стадии сокращения. Похоже, вы знаете, что делаете, поэтому я не буду подробно описывать это, но это, безусловно, важно для понимания здесь. Вы можете погуглить его и получить описания, если это необходимо. Это имеет отношение к обсуждению ниже, но я не собираюсь называть его актуальность в каждой ситуации.
Хорошо, теперь ваши вопросы:
Я хотел бы проверить мое понимание изменчивости. Volatile предотвращает неправильную оптимизацию компиляторами моего кода и обещает, что загрузка / сохранение завершаются через кеш, а не только через регистры (пожалуйста, исправьте меня, если ошиблись). Для сокращения, если суммы частичного сокращения все еще сохраняются в регистрах, почему это проблема?
Что касается определения volatile Я хотел бы направить вас к руководству по программированию CUDA. Я видел сводные описания, относящиеся к этому как к предотвращению оптимизации регистра или к переупорядочению загрузок и хранилищ. Я предпочитаю первое и буду использовать это как рабочее определение.
Основная идея заключается в том, что volatile заставляет любую ссылку (чтение или запись) на эту переменную фактически перейти в подсистему памяти. Под этим я подразумеваю, что он будет выполнять чтение или запись и не будет пытаться использовать значение, ранее загруженное в регистр. Без этого спецификатора компилятор может загрузить значение один раз (например) из фактической ячейки памяти, а затем сохранить это значение (и любые его обновления) в регистре до тех пор, пока это будет сочтено целесообразным. Компиляторы делают это с прицелом на производительность. (Кроме того, обратите внимание, что здесь вы использовали слово "кеш". Я бы избегал такого использования здесь. Общая память не имеет кеша, вставленного между ней и механизмом загрузки / сохранения процессора.)
Без volatile в этом типе синхронного кодирования мы столкнемся с проблемой, если позволим компилятору "оптимизировать" (т.е. поддерживать) промежуточные значения в регистрах. Это в первую очередь происходит из-за связи между потоками. Чтобы понять почему, давайте посмотрим на последние 2 шага вашего окончательного сокращения:
d[tid] += d[tid + 2];
d[tid] += d[tid + 1];
Давайте рассмотрим только темы, чьи tid значения 0-1. На втором-последнем шаге нить 0 подхватит d[2] значение и добавить его к d[0] значение, а поток 1 будет забрать d[3] значение и добавить его к d[1] значение. На данный момент, если мы не используем volatile компилятор не обязан писать d[1] значение, накопленное потоком 1, возвращается в общую память. Это разрешено вести в реестре. Итак d[1] значение, которое видно в разделяемой памяти, не является "актуальным".
Теперь давайте перейдем к последнему шагу. На этом этапе поток 0 читает d[1] значение из общей памяти и добавляет его к d[0] значение. Но без volatile На предыдущем этапе мы увидели, что содержимое общей памяти d[1] больше не точны. OTOH, если мы используем volatile затем произойдет запись в общую память на предыдущем шаге, а на последнем шаге поток 0 выберет правильное значение при чтении d[1], Нить CUDA - это отдельная модель. Под этим я подразумеваю, что один поток не может напрямую обращаться к значениям, содержащимся в регистрах, принадлежащих другому потоку. Таким образом, межпотоковая связь на уровне деформации обычно осуществляется либо через общую память, либо с помощью операций деформации деформации.
__syncthreads() имеет аналогичное поведение: он заставляет все оптимизированные регистром значения, подобные этим, записываться в память, чтобы они были "видны" другим потокам в блоке. Поэтому более сложной оптимизацией было бы переключиться только на volatile квалифицированный указатель, когда редукция переключается с контура __syncthreads() основанное сокращение к окончательному основному синхронному сокращению. Вы можете увидеть пример на слайдах с учебниками, на которые я ссылался в начале этого ответа.
Кроме того, в CUDA 9 такого рода синхронное программирование (более официально) не рекомендуется, вместо этого следует использовать совместные группы.
Сумма редукции даст неверные результаты без условий (n > 64) и (n > 32).
Эти условия в основном используются, потому что код разработан так, чтобы быть "правильным" для любой конфигурации блока, имеющей степень степени 2. Если мы предположим, что размер блока (количество потоков в блоке) равен степени 2 и больше 64, он должен быть, например, 128 или больше. Ваш n переменная начинается с размера блока, но затем умножается на 2:
n >>= 1;
Поэтому, если мы хотим убедиться в правильности этой строки кода:
d[tid] += d[tid + 32];
тогда мы должны применять эту операцию только тогда, когда размер блока потока равен 64 (как минимум), что все равно, что сказать, что n больше 64:
if (n > 64) d[tid] += d[tid + 32];
в отношении этого вопроса утверждается, что размещенный код ведет себя иначе, если if (n > 64) включен или нет. Причина этого заключается в том, что опубликованный код включает цикл, который пересчитывает количество потоков и количество блоков по мере сокращения:
int s = blocks;
while(s > 1) {
threads = 0;
blocks = 0;
getNumBlocksAndThreads (whichKernel, s, MAX_BLOCKS, MAX_THREADS, blocks, threads);
Этот цикл в конечном итоге приводит к размеру блока меньше 128, что означает пропуск условий if, приводящих к поломке. (просто распечатайте threads переменная, во время этого цикла).
в соответствии с этим:
У меня проблемы с поиском решения проблемы, потому что я не могу использовать функции печати, как в ЦП.
Я не уверен, в чем проблема. printf должен работать изнутри кода ядра.
Совместно используемые переменные не могут иметь инициализацию как часть их объявления согласно этому ответу. Поэтому, если n < 64, мы добавляем к сумме некоторые случайные данные массива совместно используемой памяти, что в случае ошибки.