CUDA: уменьшить алгоритм

Есть несколько вещей, на которые я бы указал в вашем коде.

1. Как общее правило / шаблон, я всегда рекомендую использовать правильная проверка ошибок cuda и запустить свой код с cuda-memcheck, в любое время у вас возникли проблемы с кодом CUDA. Тем не менее, эти методы не очень помогут с ошибкой сегмента, хотя они могут помочь позже (см. Ниже).

2. Фактическая ошибка сегмента происходит в этой строке:

   printf("sum is element is: %.f", d_out[0]);
   

   вы нарушили основное правило программирования CUDA: указатели на хосты не должны быть разыменованы в коде устройства, а указатели на устройства не должны разыменовываться в коде хоста. Это последнее условие применяется здесь. d_out указатель устройства (выделяется через cudaMalloc). Такие указатели есть нет смысла если вы попытаетесь разыменовать их в коде хоста, это приведет к ошибке сегмента.

   Решение состоит в том, чтобы скопировать данные обратно на хост до их распечатки:

   float result;
   cudaMemcpy(&result, d_out, sizeof(float), cudaMemcpyDeviceToHost);
   printf("sum is element is: %.f", result);
   

3. С помощью cudaMalloc в цикле, на той же переменной, не делая никаких cudaFree операции, не является хорошей практикой и может привести к ошибкам нехватки памяти в длительных циклах, а также может привести к программам с утечками памяти, если такая конструкция используется в более крупной программе:

   do
   {
   ...
   
   cudaMalloc(&d_intermediate, sizeof(float)*num_blocks);
   }
   while...
   

   в этом случае я думаю, что лучший подход и тривиальное решение было бы cudaFree d_intermediate прямо перед перераспределением:

   do
   {
   ...
   cudaFree(d_intermediate);
   cudaMalloc(&d_intermediate, sizeof(float)*num_blocks);
   }
   while...
   

4. Это может не делать то, что вы думаете:

   const int value = 5;
   cudaMemset(d_in, value, sizeof(float)*size);
   

   вероятно, вы знаете об этом, но cudaMemset, лайк memset, работает на количество байтов. Итак, вы заполняете d_in массив со значением, соответствующим 0x05050505 (и я понятия не имею, что соответствует этой битовой комбинации при интерпретации как float количество). Поскольку вы ссылаетесь на фиктивные ценности, вы, возможно, уже знаете об этом. Но это распространенная ошибка (например, если вы на самом деле пытались инициализировать массив со значением 5 в каждом float местоположение), поэтому я подумал, что укажу на это.

У вашего кода есть и другие проблемы (которые вы обнаружите, если внесете вышеуказанные исправления, а затем запустите свой код с cuda-memcheck). Чтобы узнать, как делать хорошие параллельные сокращения, я бы рекомендовал изучить параллельное сокращение CUDA. образец кода а также презентация. Параллельное сокращение глобальной памяти не рекомендуется по соображениям производительности.

Для полноты вот некоторые дополнительные проблемы, которые я нашел:

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

2. Ваша последняя запись в глобальную память в ядре также должна быть обусловлена ​​внутренним положением чтения. В противном случае ваша стратегия всегда запускать дополнительный блок позволила бы читать из этой строки за пределы (cuda-memcheck покажет это).

3. Логика сокращения в вашем цикле в reduce функция, как правило, испорчена и требует доработки несколькими способами.

Я не говорю, что этот код не содержит дефектов, но, похоже, он работает для данного тестового примера и дает правильный ответ (1024):

#include <stdio.h>

/* -------- KERNEL -------- */
__global__ void reduce_kernel(float * d_out, float * d_in, const int size)
{
// position and threadId
int pos = blockIdx.x * blockDim.x + threadIdx.x;
int tid = threadIdx.x;

// do reduction in global memory
for (unsigned int s = blockDim.x / 2; s>0; s>>=1)
{
if (tid < s)
{
if (pos+s < size) // Handling out of bounds
{
d_in[pos] = d_in[pos] + d_in[pos+s];
}
}
__syncthreads();
}

// only thread 0 writes result, as thread
if ((tid==0) && (pos < size))
{
d_out[blockIdx.x] = d_in[pos];
}
}

/* -------- KERNEL WRAPPER -------- */
void reduce(float * d_out, float * d_in, int size, int num_threads)
{
// setting up blocks and intermediate result holder
int num_blocks = ((size) / num_threads) + 1;
float * d_intermediate;
cudaMalloc(&d_intermediate, sizeof(float)*num_blocks);
cudaMemset(d_intermediate, 0, sizeof(float)*num_blocks);
int prev_num_blocks;
// recursively solving, will run approximately log base num_threads times.
do
{
reduce_kernel<<<num_blocks, num_threads>>>(d_intermediate, d_in, size);

// updating input to intermediate
cudaMemcpy(d_in, d_intermediate, sizeof(float)*num_blocks, cudaMemcpyDeviceToDevice);

// Updating num_blocks to reflect how many blocks we now want to compute on
prev_num_blocks = num_blocks;
num_blocks = num_blocks / num_threads + 1;

// updating intermediate
cudaFree(d_intermediate);
cudaMalloc(&d_intermediate, sizeof(float)*num_blocks);
size = num_blocks*num_threads;
}
while(num_blocks > num_threads); // if it is too small, compute rest.

// computing rest
reduce_kernel<<<1, prev_num_blocks>>>(d_out, d_in, prev_num_blocks);

}

/* -------- MAIN -------- */
int main(int argc, char **argv)
{
// Setting num_threads
int num_threads = 512;
// Making non-bogus data and setting it on the GPU
const int size = 1024;
const int size_out = 1;
float * d_in;
float * d_out;
cudaMalloc(&d_in, sizeof(float)*size);
cudaMalloc((void**)&d_out, sizeof(float)*size_out);
//const int value = 5;
//cudaMemset(d_in, value, sizeof(float)*size);
float * h_in = (float *)malloc(size*sizeof(float));
for (int i = 0; i <  size; i++) h_in[i] = 1.0f;
cudaMemcpy(d_in, h_in, sizeof(float)*size, cudaMemcpyHostToDevice);

// Running kernel wrapper
reduce(d_out, d_in, size, num_threads);
float result;
cudaMemcpy(&result, d_out, sizeof(float), cudaMemcpyDeviceToHost);
printf("sum is element is: %.f\n", result);
}