两个向量的 cuda 点积不适用于 N>=369

Question

我已经实现了一个矢量点积，如下所示。 它使用 CUDA 7.5 编译，带有 compute_20,sm_20 和 const int THREADS_PER_BLOCK = 16;。

浮点数和双精度数都会发生这种情况。

它在n=368 上有效，但除此之外，结果不正确。我想知道问题出在我的实现代码还是我正在使用的值（请参阅第二个代码，初始化），例如可能是 n=368 之外的添加引入了浮点错误（这可能很奇怪，因为浮点数和双精度数都发生了相同的错误）。

int divUp(int total, int grain) { return (total+grain-1)/grain; }

__device__ __forceinline__ double atomicAdd(double* address, double val)
{
    unsigned long long int* address_as_ull = (unsigned long long int*)address;
    unsigned long long int old = *address_as_ull, assumed;
    do
    {
        assumed = old;
        old = atomicCAS(address_as_ull, assumed, __double_as_longlong(val+__longlong_as_double(assumed)));
    }
    while(assumed!=old);
    return __longlong_as_double(old);
}

__device__ __forceinline__ float atomicAdd(float* address, float val)
{
    unsigned int *ptr = (unsigned int *)address;
    unsigned int old, newint, ret = *ptr;
    do {
        old = ret;
        newint = __float_as_int(__int_as_float(old)+val);
    } while((ret = atomicCAS(ptr, old, newint)) != old);

    return __int_as_float(ret);
}

template<typename T>
__global__ void vecdotk(const T* a, const T* b, const int n, T* c)
{
    __shared__ T temp[THREADS_PER_BLOCK];
    int x = threadIdx.x+blockIdx.x*blockDim.x;
    if(x==0) c[0] = 0.0;
    if(x<n) {temp[threadIdx.x] = a[x]*b[x];
    }
    else temp[threadIdx.x] = 0.0;
    __syncthreads();

    if(0==threadIdx.x)
    {
        T sum = 0.0;
        for(int j=0; j<THREADS_PER_BLOCK; ++j)
        {
            sum += temp[j];
        }
        atomicAdd(c, sum);
    }
}

template<typename T>
void dot(const T* a, const T* b, const int n, T* c)
{
    dim3 block(THREADS_PER_BLOCK);
    dim3 grid(divUp(n, block.x), 1);
    vecdotk<T><<<grid, block>>>(a, b, n, c);
    cudaSafeCall(cudaGetLastError());
};

我使用以下两个宿主向量来填充输入设备数组（目前我没有展示它们，因为它们是更大库的一部分）。基本上我想计算平方系列的总和，即

// fill host vectors a and b
for(int i=0; i<n; ++i)
{
    h_vec_a[i] = i+1;//__mat_rand();
    h_vec_b[i] = i+1;//__mat_rand();
}

Answer 1

这行不通：

if(x==0) c[0] = 0.0; 

不保证（在 CUDA 中）线程 0 首先运行，或者该行将在其他线程到达代码中的任何点之前运行。你需要在启动这个内核之前初始化c[0]。否则，某些线程可能会对 c 进行原子添加，然后，稍后某个时间，线程 0 可能会将 c[0] 初始化为零。

另外，CUDA 已经提供了float 版本的 atomicAdd，您没有理由提供自己的。而且，运行 16 个线程的线程块不会给你带来良好的性能（我建议只使用 CUBLAS 点积函数。）修复c[0]（删除那行代码，并在内核之前初始化c[0]）你的代码对我来说运行正确：

$ cat t372.cu
#include <stdio.h>

const int n = 2048;
#ifdef USE_DOUBLE
typedef double mytype;
#else
typedef float mytype;
#endif
const int THREADS_PER_BLOCK = 16;

int divUp(int total, int grain) { return (total+grain-1)/grain; }
#if 0
__device__ __forceinline__ double atomicAdd(double* address, double val)
{
    unsigned long long int* address_as_ull = (unsigned long long int*)address;
    unsigned long long int old = *address_as_ull, assumed;
    do
    {
        assumed = old;
        old = atomicCAS(address_as_ull, assumed, __double_as_longlong(val+__longlong_as_double(assumed)));
    }
    while(assumed!=old);
    return __longlong_as_double(old);
}

__device__ __forceinline__ float atomicAdd(float* address, float val)
{
    unsigned int *ptr = (unsigned int *)address;
    unsigned int old, newint, ret = *ptr;
    do {
        old = ret;
        newint = __float_as_int(__int_as_float(old)+val);
    } while((ret = atomicCAS(ptr, old, newint)) != old);

    return __int_as_float(ret);
}
#endif
template<typename T>
__global__ void vecdotk(const T* a, const T* b, const int n, T* c)
{
    __shared__ T temp[THREADS_PER_BLOCK];
    int x = threadIdx.x+blockIdx.x*blockDim.x;
    //if(x==0) c[0] = 0.0;
    if(x<n) {temp[threadIdx.x] = a[x]*b[x];
    }
    else temp[threadIdx.x] = 0.0;
    __syncthreads();

    if(0==threadIdx.x)
    {
        T sum = 0.0;
        for(int j=0; j<THREADS_PER_BLOCK; ++j)
        {
            sum += temp[j];
        }
        atomicAdd(c, sum);
    }
}

template<typename T>
cudaError_t dot(const T* a, const T* b, const int n, T* c)
{
    dim3 block(THREADS_PER_BLOCK);
    dim3 grid(divUp(n, block.x), 1);
    vecdotk<T><<<grid, block>>>(a, b, n, c);
    cudaDeviceSynchronize();
    return cudaGetLastError();
};

int main(){

  mytype *h_vec_a, *h_vec_b, *d_vec_a, *d_vec_b, *h_c, *d_c;
  int bs = n*sizeof(mytype);
  h_vec_a = (mytype *)malloc(bs);
  h_vec_b = (mytype *)malloc(bs);
  h_c = (mytype *)malloc(sizeof(mytype));
  cudaMalloc(&d_vec_b, bs);
  cudaMalloc(&d_vec_a, bs);
  cudaMalloc(&d_c, sizeof(mytype));
// fill host vectors a and b
  for(int i=0; i<n; ++i)
  {
    h_vec_a[i] = i+1;//__mat_rand();
    h_vec_b[i] = i+1;//__mat_rand();
  }
  h_c[0] = 0;
  cudaMemcpy(d_vec_a, h_vec_a, bs, cudaMemcpyHostToDevice);
  cudaMemcpy(d_vec_b, h_vec_b, bs, cudaMemcpyHostToDevice);
  cudaMemcpy(d_c, h_c, sizeof(mytype), cudaMemcpyHostToDevice);
  dot(d_vec_a, d_vec_b, n, d_c);
  cudaMemcpy(h_c, d_c, sizeof(mytype), cudaMemcpyDeviceToHost);
  mytype test_val = 0;
  for (int i=0; i < n; i++)
    test_val += h_vec_a[i] * h_vec_b[i];
  printf("GPU result: %f, CPU result: %f\n", h_c[0], test_val);

}
$ nvcc -arch=sm_20 -o t372 t372.cu
nvcc warning : The 'compute_20', 'sm_20', and 'sm_21' architectures are deprecated, and may be removed in a future release (Use -Wno-deprecated-gpu-targets to suppress warning).
$ cuda-memcheck ./t372
========= CUDA-MEMCHECK
GPU result: 2865411584.000000, CPU result: 2865411072.000000
========= ERROR SUMMARY: 0 errors
$

最后3位数字的差异是由于float的限制，而不是由于代码中的任何错误。例如，如果您将初始化更改为将每个向量初始化为全 1，那么在这种情况下您将获得准确的结果。

同样，出于性能原因，您的代码可能会受到许多批评。如果你想要一个快速的点产品，我建议使用CUBLAS。