CUB сокращение с использованием 2D сетки блоков

Question

Я пытаюсь сделать сумму, используя метод сокращения CUB.

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

#include <iostream>
#include <math.h>
#include <cub/block/block_reduce.cuh>
#include <cub/block/block_load.cuh>
#include <cub/block/block_store.cuh>
#include <iomanip>

#define nat 1024
#define BLOCK_SIZE 32
#define GRID_SIZE 32

struct frame
{
   int  natm;
   char  title[100];
   float conf[nat][3];
};

using namespace std;
using namespace cub;

__global__
void add(frame* s, float L, float rc, float* blocksum)
{
int i = blockDim.x*blockIdx.x + threadIdx.x;
int j = blockDim.y*blockIdx.y + threadIdx.y;

float E=0.0, rij, dx, dy, dz;

// Your calculations first so that each thread holds its result
  dx = fabs(s->conf[j][0] - s->conf[i][0]);
  dy = fabs(s->conf[j][1] - s->conf[i][1]);
  dz = fabs(s->conf[j][2] - s->conf[i][2]);
  dx = dx - round(dx/L)*L;
  dy = dy - round(dy/L)*L;
  dz = dz - round(dz/L)*L;

   rij = sqrt(dx*dx + dy*dy + dz*dz);

  if ((rij <= rc) && (rij > 0.0))
    {E =  (4*((1/pow(rij,12))-(1/pow(rij,6))));}

//  E = 1.0;
__syncthreads();
// Block wise reduction so that one thread in each block holds sum of thread results

typedef cub::BlockReduce<float, BLOCK_SIZE, BLOCK_REDUCE_RAKING, BLOCK_SIZE> BlockReduce;

__shared__ typename BlockReduce::TempStorage temp_storage;

float aggregate = BlockReduce(temp_storage).Sum(E);

if (threadIdx.x == 0 && threadIdx.y == 0)
    blocksum[blockIdx.x*blockDim.y + blockIdx.y] = aggregate;

}

int main(void)
{
  frame  * state = (frame*)malloc(sizeof(frame));

  float *blocksum = (float*)malloc(GRID_SIZE*GRID_SIZE*sizeof(float));

  state->natm = nat; //inicializando o numero de atomos;

  char name[] = "estado1";
  strcpy(state->title,name);

  for (int i = 0; i < nat; i++) {
    state->conf[i][0] = i;
    state->conf[i][1] = i;
    state->conf[i][2] = i;
  }

  frame * d_state;
  float *d_blocksum;

  cudaMalloc((void**)&d_state, sizeof(frame));

  cudaMalloc((void**)&d_blocksum, ((GRID_SIZE*GRID_SIZE)*sizeof(float)));

  cudaMemcpy(d_state, state, sizeof(frame),cudaMemcpyHostToDevice);


  dim3 dimBlock(BLOCK_SIZE,BLOCK_SIZE);
  dim3 gridBlock(GRID_SIZE,GRID_SIZE);

  add<<<gridBlock,dimBlock>>>(d_state, 3000, 15, d_blocksum);

  cudaError_t status =  cudaMemcpy(blocksum, d_blocksum, ((GRID_SIZE*GRID_SIZE)*sizeof(float)),cudaMemcpyDeviceToHost);

  float Etotal = 0.0;
  for (int k = 0; k < GRID_SIZE*GRID_SIZE; k++){
       Etotal += blocksum[k];
  }
 cout << endl << "energy: " << Etotal << endl;

  if (cudaSuccess != status)
  {
    cout << cudaGetErrorString(status) << endl;
  }

 // Free memory
  cudaFree(d_state);
  cudaFree(d_blocksum);

  return cudaThreadExit();
}

Что происходит, что если значение GRID_SIZE такой же какBLOCK_SIZE, как написано выше. Расчет правильный. Но если я изменю значение GRID_SIZE, результат идет не так, как надо. Что заставляет меня думать, что ошибка в этом коде:

blocksum[blockIdx.x*blockDim.y + blockIdx.y] = aggregate;

Идея здесь состоит в том, чтобы вернуть одномерный массив, который содержит сумму каждого блока.

Я не собираюсь менять BLOCK_SIZE значение, но значение GRID_SIZE В зависимости от системы, которую я смотрю, я намерен использовать значения больше 32 (всегда кратно этому).

Я искал пример, который использует 2D-сетку с CUB, но не нашел.

Я действительно новичок в программе CUDA, возможно, я делаю ошибку.

редактировать: я поставил полный код. Для сравнения, когда я вычисляю эти точные значения для последовательной программы, это дает мне энергию: -297,121

Answer 1

Вероятно, главная проблема в том, что ваша выходная индексация не верна. Вот сокращенная версия вашего кода, демонстрирующая правильные результаты для произвольного GRID_SIZE:

$ cat t1360.cu
#include <stdio.h>
#include <cub/cub.cuh>
#define BLOCK_SIZE 32
#define GRID_SIZE 25
__global__
void add(float* blocksum)
{
   float E = 1.0;
  // Block wise reduction so that one thread in each block holds sum of thread results
    typedef cub::BlockReduce<float, BLOCK_SIZE, cub::BLOCK_REDUCE_RAKING, BLOCK_SIZE> BlockReduce;

    __shared__ typename BlockReduce::TempStorage temp_storage;
    float aggregate = BlockReduce(temp_storage).Sum(E);
    __syncthreads();
    if (threadIdx.x == 0 && threadIdx.y == 0)
        blocksum[blockIdx.y*gridDim.x + blockIdx.x] = aggregate;
}

int main(){

  float *d_result, *h_result;
  h_result = (float *)malloc(GRID_SIZE*GRID_SIZE*sizeof(float));
  cudaMalloc(&d_result, GRID_SIZE*GRID_SIZE*sizeof(float));
  dim3 grid  = dim3(GRID_SIZE,GRID_SIZE);
  dim3 block = dim3(BLOCK_SIZE, BLOCK_SIZE);
  add<<<grid, block>>>(d_result);
  cudaMemcpy(h_result, d_result, GRID_SIZE*GRID_SIZE*sizeof(float), cudaMemcpyDeviceToHost);
  cudaError_t err = cudaGetLastError();
  if (err != cudaSuccess) {printf("cuda error: %s\n", cudaGetErrorString(err)); return -1;}
  float result = 0;
  for (int i = 0; i < GRID_SIZE*GRID_SIZE; i++) result += h_result[i];
  if (result != (float)(GRID_SIZE*GRID_SIZE*BLOCK_SIZE*BLOCK_SIZE)) printf("mismatch, should be: %f, was: %f\n", (float)(GRID_SIZE*GRID_SIZE*BLOCK_SIZE*BLOCK_SIZE), result);
  else printf("Success\n");
  return 0;
}

$ nvcc -o t1360 t1360.cu
$ ./t1360
Success
$

Важное изменение, которое я внес в код вашего ядра, заключалось в индексировании выходных данных:

blocksum[blockIdx.y*gridDim.x + blockIdx.x] = aggregate;

Мы хотим смоделировать 2D-индекс в массив, который имеет ширину и высоту GRID_SIZE состоящий из одного float количество на балл. Поэтому ширина этого массива определяется как gridDim.x (не blockDim). gridDim Переменная дает размеры сетки в терминах блоков - и это точно соответствует тому, как настроен наш массив результатов.

Ваш опубликованный код потерпит неудачу, если GRID_SIZE а также BLOCK_SIZE разные (например, если GRID_SIZE были меньше чем BLOCK_SIZE, cuda-memcheck покажет незаконный доступ, и если GRID_SIZE больше чем BLOCK_SIZE тогда эта ошибка индексации приведет к тому, что блоки перезаписывают значения друг друга в выходном массиве) из-за этой путаницы между blockDim а также gridDim,

Также обратите внимание, что float Операции обычно имеют только около 5 десятичных цифр точности. Столь небольшие различия в 5-м или 6-м десятичном знаке могут быть связаны с разницей в порядке операций при выполнении арифметики с плавающей запятой. Вы можете доказать это себе, переключившись на double арифметика.