Реализация Max Reduce в Cuda

Ваш основной вывод в вашем удаленном ответе был правильным: опубликованное вами ядро ​​​​не понимает того факта, что в конце выполнения этого ядра вы выполнили значительную часть общего сокращения, но результаты не совсем полные. Результаты каждого блока должны быть объединены (каким-то образом). Как указано в комментариях, есть еще несколько проблем с вашим кодом. Давайте посмотрим на его модифицированную версию:

__device__ float atomicMaxf(float* address, float val)
{
    int *address_as_int =(int*)address;
    int old = *address_as_int, assumed;
    while (val > __int_as_float(old)) {
        assumed = old;
        old = atomicCAS(address_as_int, assumed,
                        __float_as_int(val));
        }
    return __int_as_float(old);
}


__global__ void max_reduce(const float* const d_array, float* d_max, 
                                              const size_t elements)
{
    extern __shared__ float shared[];

    int tid = threadIdx.x;
    int gid = (blockDim.x * blockIdx.x) + tid;
    shared[tid] = -FLOAT_MAX;  // 1

    if (gid < elements)
        shared[tid] = d_array[gid];
    __syncthreads();

    for (unsigned int s=blockDim.x/2; s>0; s>>=1) 
    {
        if (tid < s && gid < elements)
            shared[tid] = max(shared[tid], shared[tid + s]);  // 2
        __syncthreads();
    }
    // what to do now?
    // option 1: save block result and launch another kernel
    if (tid == 0)        
        d_max[blockIdx.x] = shared[tid]; // 3
    // option 2: use atomics
    if (tid == 0)
      atomicMaxf(d_max, shared[0]);
}

1. Как указал Паван, вам необходимо инициализировать массив общей памяти. Последний запущенный блок может не быть «полным», если gridDim.x*blockDim.x больше elements.
2. Обратите внимание, что в этой строке, даже несмотря на то, что мы проверяем, что рабочий поток (gid) меньше elements, когда мы добавляем s к gid для индексации в разделяемую память, мы все еще можем индексировать за пределами легитимного значения, скопированные в разделяемую память, в последнем блоке. Поэтому нам нужна инициализация разделяемой памяти, указанная в примечании 1.
3. Как вы уже обнаружили, ваша последняя строка была неверной. Каждый блок выдает свой результат, и мы должны их как-то комбинировать. Один из методов, который вы можете рассмотреть, если количество запущенных блоков невелико (подробнее об этом позже), заключается в использовании атомарность. Обычно мы отговариваем людей от использования атомарных вычислений, поскольку они «затратны» с точки зрения времени выполнения. Однако другой вариант, с которым мы сталкиваемся, — это сохранение результата блока в глобальной памяти, завершение работы ядра, а затем, возможно, запуск другого ядра для объединения результатов отдельных блоков. Если я изначально запустил большое количество блоков (скажем, более 1024), то, если я буду следовать этой методологии, я могу в конечном итоге запустить два дополнительных ядра. Таким образом, рассмотрение атомарности. Как указано, встроенной функции atomicMax для чисел с плавающей запятой нет, но, как указано в документации, вы можете использовать atomicCAS для создания любой произвольной атомарной функции, и я привел пример в atomicMaxf, который обеспечивает атомарный максимум для float.

Но лучше ли запускать 1024 или более атомарных функций (по одной на блок)? Возможно нет.

При запуске ядра блоков потоков нам действительно нужно запустить достаточное количество блоков потоков, чтобы компьютер был занят. Как правило, мы хотим, чтобы на один SM работало не менее 4-8 варпов, и, вероятно, неплохо было бы сделать несколько больше. Но с точки зрения использования машины нет особой выгоды от первоначального запуска тысяч блоков потоков. Если мы выберем число, например, 8 блоков потоков на SM, а у нас самое большее, скажем, 14-16 SM в нашем графическом процессоре, это даст нам относительно небольшое число 8*14 = 112 блоков потоков. Давайте выберем 128 (8*16) для красивого круглого числа. В этом нет ничего волшебного, этого достаточно, чтобы GPU был загружен. Если мы заставим каждый из этих 128 блоков выполнять дополнительную работу для решения всей проблемы, мы сможем использовать атомарность без (возможно) больших штрафов за это и избежать множественного ядра. запускает. Итак, как это будет выглядеть?:

__device__ float atomicMaxf(float* address, float val)
{
    int *address_as_int =(int*)address;
    int old = *address_as_int, assumed;
    while (val > __int_as_float(old)) {
        assumed = old;
        old = atomicCAS(address_as_int, assumed,
                        __float_as_int(val));
        }
    return __int_as_float(old);
}


__global__ void max_reduce(const float* const d_array, float* d_max, 
                                              const size_t elements)
{
    extern __shared__ float shared[];

    int tid = threadIdx.x;
    int gid = (blockDim.x * blockIdx.x) + tid;
    shared[tid] = -FLOAT_MAX; 

    while (gid < elements) {
        shared[tid] = max(shared[tid], d_array[gid]);
        gid += gridDim.x*blockDim.x;
        }
    __syncthreads();
    gid = (blockDim.x * blockIdx.x) + tid;  // 1
    for (unsigned int s=blockDim.x/2; s>0; s>>=1) 
    {
        if (tid < s && gid < elements)
            shared[tid] = max(shared[tid], shared[tid + s]);
        __syncthreads();
    }

    if (tid == 0)
      atomicMaxf(d_max, shared[0]);
}

С этим модифицированным ядром при создании запуска ядра мы не решаем, сколько блоков потоков запускать, исходя из общего размера данных (elements). Вместо этого мы запускаем фиксированное количество блоков (скажем, 128, вы можете изменить это число, чтобы узнать, что работает быстрее) и позволяем каждому блоку потока (и, следовательно, всей сетке) циклически проходить по памяти, вычисляя частичные максимальные операции над каждым элементом в Общая память. Затем в строке, отмеченной комментарием 1, мы должны переустановить переменную gid в исходное значение. Это на самом деле не нужно, и код цикла сокращения блоков можно еще больше упростить, если мы гарантируем, что размер сетки (gridDim.x*blockDim.x) меньше, чем elements, что несложно сделать при запуске ядра.

Обратите внимание, что при использовании этого атомарного метода необходимо инициализировать результат (в данном случае *d_max) соответствующим значением, например -FLOAT_MAX.

Опять же, обычно мы уводим людей от атомарного использования, но в этом случае стоит подумать, тщательно ли мы управляем им, и это позволяет нам сэкономить накладные расходы на дополнительный запуск ядра.

Для анализа на уровне ниндзя того, как выполнять быстрые параллельные сокращения, взгляните на отличный технический документ Марка Харриса, который доступен с соответствующим Пример CUDA.

Вот тот, который кажется наивным, но это не так. Это не будет распространяться на другие функции, такие как sum(), но отлично работает для min() и max().

__device__ const float float_min = -3.402e+38;

__global__ void maxKernel(float* d_data)
{ 
    // compute max over all threads, store max in d_data[0]
    int i = threadIdx.x;
    __shared__ float max_value;

    if (i == 0) max_value = float_min;
    float v = d_data[i];
    __syncthreads();

    while (max_value < v) max_value = v;

    __syncthreads();
    if (i == 0) d_data[0] = max_value;
}

Да, правильно, только один раз синхронизируется после инициализации и один раз перед записью результата. К черту условия гонки! Полный вперед!

Прежде чем вы скажете мне, что это не сработает, пожалуйста, сначала попробуйте. Я тщательно тестировал, и он каждый раз работает на множестве произвольных размеров ядра. Оказывается, в данном случае условие гонки не имеет значения, потому что его разрешает цикл while.

Он работает значительно быстрее, чем обычное сокращение. Еще одним сюрпризом является то, что среднее количество проходов для ядра размером 32 равно 4. Да, это (log(n)-1), что кажется нелогичным. Это потому, что состояние гонки дает шанс на удачу. Этот бонус предоставляется в дополнение к устранению накладных расходов, связанных с обычным сокращением.

При большем n невозможно избежать по крайней мере одной итерации на деформацию, но эта итерация включает только одну операцию сравнения, которая обычно сразу ложна по всей деформации, когда max_value находится на верхнем конце распределения. Вы можете изменить его, чтобы использовать несколько SM, но это значительно увеличит общую рабочую нагрузку и увеличит стоимость связи, поэтому вряд ли поможет.

Для краткости я опустил размер и выходные аргументы. Размер — это просто количество потоков (которое может быть 137 или сколько угодно). Вывод возвращается в d_data[0].

Я загрузил рабочий файл здесь: https://github.com/kenseehart/YAMR.