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, возможно, я делаю ошибку.

edit : Я поставил полный код.Для сравнения, когда я вычисляю эти точные значения для последовательной программы, это дает мне энергию: -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;

Мыхотите смоделировать двумерный индекс в массив, который имеет ширину и высоту 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 арифметику.