0

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

#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 reutrns 6.73066e+07 而 MATLAB 返回 9.3823e+07。对于 N=2048*64 CUDA 给出 3.28033e+08 而 MATLAB给出 7.5059e+08)。我倾向于相信差异源于float在我的 C 代码中的使用，但如果我用编译器替换它，则会double抱怨atomicAdd不支持双参数。我应该如何解决这个问题？

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

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

4

1 回答 1

2

从内核中注释掉这一行：

//   *c = 0.00;

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

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

而且我相信您会或多或少地获得与 matlab 匹配的结果。

事实上，你的内核中的这条线不受任何同步的保护，这让你很困惑。每个块的每个线程，无论何时执行，都会c在您编写它时归零。

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

编辑：在评论中回答以下问题：

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

于 2013-04-04T22:34:18.947 回答