CUDA - почему параллельное сокращение, основанное на деформации, замедляется?

У меня возникла идея о параллельном сокращении, основанном на warp, поскольку все потоки warp синхронизируются по определению.

Итак, идея заключалась в том, что входные данные могут быть уменьшены в 64 раз (каждый поток уменьшает два элемента) без какой-либо необходимости синхронизации.

То же, что и оригинальная реализация Марк Харрис, сокращение применяется на уровне блока, а данные - в общей памяти. http://gpgpu.org/static/sc2007/SC07_CUDA_5_Optimization_Harris.pdf

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

Сам алгоритм отлично работает. Протестировано с полным массивом, чтобы проверить "подсчет".

Тело функции реализаций:

/**
 * Performs a parallel reduction with operator add 
 * on the given array and writes the result with the thread 0
 * to the given target value
 *
 * @param inValues T* Input float array, length must be a multiple of 2 and equal to blockDim.x
 * @param targetValue float 
 */
__device__ void reductionAddBlockThread_f(float* inValues,
    float &outTargetVar)
{
    // code of the below functions
}

1. Реализация его версии:

if (blockDim.x >= 1024 && threadIdx.x < 512)
    inValues[threadIdx.x] += inValues[threadIdx.x + 512];
__syncthreads();
if (blockDim.x >= 512 && threadIdx.x < 256)
    inValues[threadIdx.x] += inValues[threadIdx.x + 256];
__syncthreads();
if (blockDim.x >= 256 && threadIdx.x < 128)
    inValues[threadIdx.x] += inValues[threadIdx.x + 128];
__syncthreads();
if (blockDim.x >= 128 && threadIdx.x < 64)
    inValues[threadIdx.x] += inValues[threadIdx.x + 64];
__syncthreads();

//unroll last warp no sync needed
if (threadIdx.x < 32)
{
    if (blockDim.x >= 64) inValues[threadIdx.x] += inValues[threadIdx.x + 32];
    if (blockDim.x >= 32) inValues[threadIdx.x] += inValues[threadIdx.x + 16];
    if (blockDim.x >= 16) inValues[threadIdx.x] += inValues[threadIdx.x + 8];
    if (blockDim.x >= 8) inValues[threadIdx.x] += inValues[threadIdx.x + 4];
    if (blockDim.x >= 4) inValues[threadIdx.x] += inValues[threadIdx.x + 2];
    if (blockDim.x >= 2) inValues[threadIdx.x] += inValues[threadIdx.x + 1];

    //set final value
    if (threadIdx.x == 0)
        outTargetVar = inValues[0];
}

Ressources:

4 используемых syncthreads 12, если утверждения используются
11 читать + добавлять + писать операции
1 окончательная операция записи
5 использование регистра

Производительность:

пять тестовых пробегов: ~ 19,54 мс

2. Основанный на Warp подход: (Те же функциональные тела, что и выше)

/*
 * Perform first warp based reduction by factor of 64
 *
 * 32 Threads per Warp -> LOG2(32) = 5
 *
 * 1024 Threads / 32 Threads per Warp = 32 warps
 * 2 elements compared per thread -> 32 * 2 = 64 elements per warp
 *
 * 1024 Threads/elements divided by 64 = 16
 * 
 * Only half the warps/threads are active
 */
if (threadIdx.x < blockDim.x >> 1)
{
    const unsigned int warpId = threadIdx.x >> 5;
    // alternative threadIdx.x & 31
    const unsigned int threadWarpId = threadIdx.x - (warpId << 5);
    const unsigned int threadWarpOffset = (warpId << 6) + threadWarpId;

    inValues[threadWarpOffset] += inValues[threadWarpOffset + 32];
    inValues[threadWarpOffset] += inValues[threadWarpOffset + 16];
    inValues[threadWarpOffset] += inValues[threadWarpOffset + 8];
    inValues[threadWarpOffset] += inValues[threadWarpOffset + 4];
    inValues[threadWarpOffset] += inValues[threadWarpOffset + 2];
    inValues[threadWarpOffset] += inValues[threadWarpOffset + 1];
}

// synchronize all warps - the local warp result is stored
// at the index of the warp equals the first thread of the warp
__syncthreads();

// use first warp to reduce the 16 warp results to the final one
if (threadIdx.x < 8)
{
    // get first element of a warp
    const unsigned int warpIdx = threadIdx.x << 6;

    if (blockDim.x >= 1024) inValues[warpIdx] += inValues[warpIdx + 512];
    if (blockDim.x >= 512) inValues[warpIdx] += inValues[warpIdx + 256];
    if (blockDim.x >= 256) inValues[warpIdx] += inValues[warpIdx + 128];
    if (blockDim.x >= 128) inValues[warpIdx] += inValues[warpIdx + 64];

    //set final value
    if (threadIdx.x == 0)
        outTargetVar = inValues[0];
}

Ressources:

1 используемый syncthread 7 если утверждения
10 читать добавить операции записи
1 окончательная операция записи
5 использование регистра

5-битные сдвиги
1 добавить 1 sub

Производительность:

пять тестовых пробегов: ~ 20,82 мс

Тестирование обоих ядер несколько раз на Geforce 8800 GT 512 МБ с 256 mb значений float. И работающее ядро с 256 потоками на блок (100% заполняемость).

версия на основе warp ~ 1,28 миллисекунды медленнее.

Если будущая карта допускает большие размеры блоков, основанный на основах подход все равно не нуждается в дополнительной инструкции синхронизации, так как max 4096, который уменьшается до 64, которые уменьшаются на окончательный warp до 1

Почему это не быстрее?, или где ошибка в идее, ядро?

Из использования ressources подход warp должен быть впереди?

Edit1: Исправлено ядро, в котором активна только половина потоков, что не приводит к ошибкам чтения, добавляет новые данные о производительности

Ответ 1

Я думаю, причина, по которой ваш код медленнее, чем мой, заключается в том, что в моем коде на первом этапе в каждом ADD задействовано пополам столько же перекосов. В вашем коде все перекосы активны для всего первого этапа. Таким образом, ваш код выполняет больше команд warp. В CUDA важно учитывать тотальные инструкции "warp", а не только количество команд, выполняемых одним варпом.

Кроме того, нет смысла использовать только половину ваших перекосов. Есть накладные расходы на запуск перекосов только для того, чтобы они оценивали две ветки и выходили.

Другая мысль заключается в том, что использование unsigned char и short может фактически стоить вам производительности. Я не уверен, но это, безусловно, не сохраняет ваши регистры, поскольку они не упакованы в одиночные 32-битные переменные.

Кроме того, в моем исходном коде я заменил blockDim.x параметром шаблона BLOCKDIM, что означает, что он использовал только 5 операторов времени выполнения (ifs на втором этапе устраняются компилятором).

Кстати, более дешевый способ вычисления вашего threadWarpId -

const int threadWarpId = threadIdx.x & 31;

Вы можете проверить эту статью для получения дополнительных идей.

EDIT: Здесь альтернативное сокращение блоков на основе основы.

template <typename T, int level>
__device__
void sumReduceWarp(volatile T *sdata, const unsigned int tid)
{
  T t = sdata[tid];
  if (level > 5) sdata[tid] = t = t + sdata[tid + 32];
  if (level > 4) sdata[tid] = t = t + sdata[tid + 16];
  if (level > 3) sdata[tid] = t = t + sdata[tid +  8];
  if (level > 2) sdata[tid] = t = t + sdata[tid +  4];
  if (level > 1) sdata[tid] = t = t + sdata[tid +  2];
  if (level > 0) sdata[tid] = t = t + sdata[tid +  1];
}

template <typename T>
__device__
void sumReduceBlock(T *output, volatile T *sdata)
{
  // sdata is a shared array of length 2 * blockDim.x

  const unsigned int warp = threadIdx.x >> 5;
  const unsigned int lane = threadIdx.x & 31;
  const unsigned int tid  = (warp << 6) + lane;

  sumReduceWarp<T, 5>(sdata, tid);
  __syncthreads();

  // lane 0 of each warp now contains the sum of two warp values
  if (lane == 0) sdata[warp] = sdata[tid];

  __syncthreads();

  if (warp == 0) {
    sumReduceWarp<T, 4>(sdata, threadIdx.x);
    if (lane == 0) *output = sdata[0];
  }
}

Это должно быть немного быстрее, потому что он использует все начальные напряжения, которые запускаются на первом этапе, и не имеет ветвления на последнем этапе за счет дополнительной ветки, общей загрузки/хранения и __syncthreads() в новый средний stage. Я не тестировал этот код. Если вы запустите его, дайте мне знать, как он работает. Если вы используете шаблон для blockDim в своем исходном коде, он может быть быстрее, но я думаю, что этот код более краткий.

Обратите внимание, что временная переменная t используется, потому что Fermi и более поздние архитектуры используют чистую архитектуру загрузки/хранения, поэтому += из общей памяти в общую память приводит к дополнительной нагрузке (поскольку указатель sdata должен быть изменчивым). Явно загружаю во временное однократно, избегая этого. На G80 это не повлияет на производительность.

Ответ 1

Ответ 2