使用原子操作的 CUDA 中的点积 - 得到错误的结果

Posted 2023-02-16

【中文标题】使用原子操作的 CUDA 中的点积 - 得到错误的结果

【英文标题】：Dot Product in CUDA using atomic operations - getting wrong results

【发布时间】：2013-04-04 21:52:10

【问题描述】：

我正在尝试在 CUDA 中实现点积，并将结果与​​ MATLAB 返回的结果进行比较。我的 CUDA 代码（基于this tutorial）如下：

#include <stdio.h>

#define N (2048 * 8)
#define THREADS_PER_BLOCK 512
#define num_t float

// The kernel - DOT PRODUCT
__global__ void dot(num_t *a, num_t *b, num_t *c) 

  __shared__ num_t temp[THREADS_PER_BLOCK];
  int index = threadIdx.x + blockIdx.x * blockDim.x;
  temp[threadIdx.x] = a[index] * b[index];
  __syncthreads(); //Synchronize!
  *c = 0.00;
  // Does it need to be tid==0 that
  // undertakes this task?
  if (0 == threadIdx.x) 
    num_t sum = 0.00;
    int i;
    for (i=0; i<THREADS_PER_BLOCK; i++)
      sum += temp[i];
    atomicAdd(c, sum);        
    //WRONG: *c += sum; This read-write operation must be atomic!
  



// Initialize the vectors:
void init_vector(num_t *x)

  int i;
  for (i=0 ; i<N ; i++)
    x[i] = 0.001 * i;
  


// MAIN
int main(void)

  num_t *a, *b, *c;
  num_t *dev_a, *dev_b, *dev_c;
  size_t size = N * sizeof(num_t);

  cudaMalloc((void**)&dev_a, size);
  cudaMalloc((void**)&dev_b, size);
  cudaMalloc((void**)&dev_c, size);

  a = (num_t*)malloc(size);
  b = (num_t*)malloc(size);
  c = (num_t*)malloc(size);

  init_vector(a);
  init_vector(b);

  cudaMemcpy(dev_a, a, size, cudaMemcpyHostToDevice);
  cudaMemcpy(dev_b, b, size, cudaMemcpyHostToDevice);

  dot<<<N/THREADS_PER_BLOCK, THREADS_PER_BLOCK>>>(dev_a, dev_b, dev_c);

  cudaMemcpy(c, dev_c, sizeof(num_t), cudaMemcpyDeviceToHost);

  printf("a = [\n");
  int i;
  for (i=0;i<10;i++)
    printf("%g\n",a[i]);
  
  printf("...\n");
  for (i=N-10;i<N;i++)
    printf("%g\n",a[i]);
  
  printf("]\n\n");
  printf("a*b = %g.\n", *c);


  free(a); free(b); free(c);

  cudaFree(dev_a);
  cudaFree(dev_b);
  cudaFree(dev_c);

然后我编译它：

/usr/local/cuda-5.0/bin/nvcc -m64  -I/usr/local/cuda-5.0/include -gencode arch=compute_20,code=sm_20 -o multi_dot_product.o -c multi_dot_product.cu
g++ -m64 -o multi_dot_product multi_dot_product.o -L/usr/local/cuda-5.0/lib64 -lcudart

有关我的 NVIDIA 显卡的信息，请访问 http://pastebin.com/8yTzXUuK。我尝试使用以下简单代码在 MATLAB 中验证结果：

N = 2048 * 8;
a = zeros(N,1);
for i=1:N
    a(i) = 0.001*(i-1);
end

dot_product = a'*a;

但是随着 N 的增加，我得到了显着不同的结果（例如，对于 N=2048*32，CUDA 返回 6.73066e+07，而 MATLAB 返回 9.3823e+07。对于 N=2048*64，CUDA 给出 3.28033e+ 08 而 MATLAB 给出 7.5059e+08)。我倾向于相信差异源于在我的 C 代码中使用 float，但如果我用 double 替换它，编译器会抱怨 atomicAdd 不支持双参数。我该如何解决这个问题？

更新：另外，对于 N 的高值（例如 2048*64），我注意到 CUDA 返回的结果在每次运行时都会发生变化。如果 N 较低（例如 2048*8），则不会发生这种情况。

同时我有一个更基本的问题：变量temp 是一个大小为THREADS_PER_BLOCK 的数组，并且在同一块中的线程之间共享。它是否也在块之间共享，或者每个块都对该变量的不同副本进行操作？我是否应该将dot 方法视为对每个块的说明？有人可以详细说明在这个例子中工作是如何分割的以及变量是如何共享的

【问题讨论】：

你能量化什么是“显着不同的结果”吗？什么是相对误差和绝对误差？

@talonmies 是的，我更新了我的问题。

每个块都在任何给定 __shared__ 变量的不同副本上运行，包括您的 temp。您可以使用here 概述的方法执行double atomicAdd。

你有这行：*c = 0.00; 不受内核代码中任何同步保护的事实让你一团糟。在调用内核之前将c 初始化为零，并从内核中删除该行。每个块的每个线程，无论何时执行，都会将 c 归零，就像您编写的那样。

【参考方案1】：

从你的内核中注释掉这一行：

//   *c = 0.00;

并将这些行添加到您的主机代码中，在内核调用之前（dev_c 的 cudaMalloc 之后）：

num_t h_c = 0.0f;
cudaMemcpy(dev_c, &h_c, sizeof(num_t), cudaMemcpyHostToDevice);

而且我相信你会得到或多或少与 matlab 匹配的结果。

您的内核中有这一行不受任何同步保护的事实让您感到困惑。每个块的每个线程，无论何时执行，都会将 c 归零，就像您编写的那样。

顺便说一句，通过使用经典的并行归约方法，我们通常可以在此操作上做得更好。基本（未优化）插图是here。如果您将该方法与共享内存的使用和最后的单个 atomicAdd（每个块一个 atomicAdd）结合起来，您将获得显着改进的实现。虽然它不是点积，但this example 结合了这些想法。

编辑：在 cmets 中回答以下问题：

内核函数是网格中的所有线程（根据定义，与内核启动相关的所有线程）执行的指令集。但是，将执行视为由线程块管理是合理的，因为线程块中的线程在很大程度上是一起执行的。然而，即使在一个线程块内，执行也不一定在所有线程之间完美同步。通常，当我们想到锁步执行时，我们会想到 warp，它是单个线程块中的一组 32 个线程。因此，由于块内的扭曲之间的执行可能会出现偏差，因此即使对于单个线程块也存在这种危险。但是，如果只有一个线程块，我们可以使用适当的同步和控制机制（如 __syncthreads() 和 (if threadIdx.x == 0) 等）消除代码中的危险。但是这些机制对于控制跨多个执行的一般情况是无用的线程块。多个线程块可以按任何顺序执行。整个网格中唯一定义的同步机制是内核启动本身。因此，为了解决您的问题，我们必须在内核启动之前将 c 归零。

【讨论】：

该死的，对！因此，变量的初始化应该以受保护/同步的方式完成（就像我对 dev_a 和 dev_b 所做的那样）。还有一个问题要确保我做对了：内核函数应该被视为要在一个块中执行的一堆指令，对吗？所以，这就是我们使用temp[threadIdx.x]的原因。

在上面的答案中回答了这个问题