Ошибка суммы префикса CUDA Paralell

Я пытаюсь реализовать трехфазное параллельное сканирование, как описано в главе 8 "Программирование массивно-параллельных процессоров" 3-го издания (есть какая-то строка кода, но только инструкции). Этот алгоритм позволяет использовать только 1 блок с максимальным количеством потоков в блоке, и он ограничен размером разделяемой памяти, поскольку все элементы должны помещаться в разделяемую память

После некоторой отладки я обнаружил проблему во время суммирования в фазе 3 (см. Код ниже), когда я использую большое количество элементов, например 8192, и более 1 потока.

Графическая концепция алгоритма заключается в следующем: введите описание изображения здесь

А ниже вы можете увидеть код ядра:

__global__ 
void efficient_Kogge_Stone_scan_kernel(float *X, float *Y, int InputSize) {
    __shared__ float XY[SECTION_SIZE];
    __shared__ float AUS[BLOCK_DIM];
    //int i = blockIdx.x * blockDim.x + threadIdx.x;

    // Keep mind: Partition the input into blockDim.x subsections: i.e. for 8 threads --> 8 subsections

    // collaborative load in a coalesced manner
    for (int j = 0; j < SECTION_SIZE; j += blockDim.x) {
        XY[threadIdx.x + j] = X[threadIdx.x + j];
    }
    __syncthreads();


    // PHASE 1: scan inner own subsection
    // At the end of this phase the last element of each subsection contains the sum of all alements in own subsection
    for (int j = 1; j < SUBSECTION_SIZE; j++) {
        XY[threadIdx.x * (SUBSECTION_SIZE)+j] += XY[threadIdx.x * (SUBSECTION_SIZE)+j - 1];
    }
    __syncthreads();


    // PHASE 2: perform iterative kogge_stone_scan of the last elements of each subsections of XY loaded first in AUS
    AUS[threadIdx.x] = XY[threadIdx.x * (SUBSECTION_SIZE)+(SUBSECTION_SIZE)-1];
    float in;
    for (unsigned int stride = 1; stride < blockDim.x; stride *= 2) {
        __syncthreads();
        if (threadIdx.x >= stride) {
            in = AUS[threadIdx.x - stride];
        }
        __syncthreads();
        if (threadIdx.x >= stride) {
            AUS[threadIdx.x] += in;
        }
    }
    __syncthreads();


    // PHASE 3: each thread adds to its elements the new value of the last element of its predecessor's section
    if (threadIdx.x > 0) {
        for (unsigned int stride = 0; stride < (SUBSECTION_SIZE); stride++) {
            XY[threadIdx.x * (SUBSECTION_SIZE)+stride] += AUS[threadIdx.x - 1];  // <-- 
        }
    }
    __syncthreads();


    // store the result into output vector
    for (int j = 0; j < SECTION_SIZE; j += blockDim.x) {
        Y[threadIdx.x + j] = XY[threadIdx.x + j];
    }
}

Если я использую один поток в блоке и 8192 элемента, он работает отлично, но если я использую более одного потока, результат будет неправильным в XY[5793] (или X[5793], когда он закончил и сохранил результаты). Он работает с 4096 элементами и одним или несколькими потоками до 1024 потоков. Если я использую int вместо чисел с плавающей точкой, он работает даже с 8192 элементами с одним или несколькими потоками.

Я попытался проверить также в MATLAB, и это сравнение результатов:

  • X [5973] = 16788115 ---- MATLAB
  • X [5973] = 16788114 ---- CPU
  • X [5973] = 16788116 ---- GPU

Как мы видим, результат процессора также отличается от MATLAB, поэтому после этих результатов я думаю, что проблема заключается в сложении с плавающей запятой, но я говорю вам, что я заполнил входной массив упорядоченными числами с плавающей запятой "x.00" (например, {1,00, 2,00, 3,00, 4,00 ..... 8192,00}).

Другая проблема связана с производительностью, код хоста всегда быстрее, чем код ядра, с этими параметрами конфигурации и этими входами, это нормально?

Если вам нужен весь исходный код, вы можете найти его здесь

Источник

2 ответа

Решение

8192 - это 2^13

сумма (1..8192) составляет около 8192^2/2: 8192*8193/2, то есть чуть больше 2^25.

Таким образом, вам нужно 26 бит для его представления (см. Примечание ниже).

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

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

Источник

Может быть проблема с первым сканированием:

XY[threadIdx.x * (SUBSECTION_SIZE)+j] += XY[threadIdx.x * (SUBSECTION_SIZE)+j - 1];

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

Попробуйте разбить этот раздел на части, сохранив значение в регистре. Пример:

int t =  XY[threadIdx.x * (SUBSECTION_SIZE)+j - 1];
 __syncthreads();
 XY[threadIdx.x * (SUBSECTION_SIZE)+j] += t;
Источник
Другие вопросы по тегам