CUDA Dot产品

Question

我正在尝试实现双精度数组的经典点积内核，并通过各种块的最终总和的原子计算。我使用atomicAdd作为编程指南的第116页所述的双精度。可能是我做错了什么。每个块中线程的部分总和计算正确，但后缀的原子操作似乎没有正常工作因为每次我用相同的数据运行我的内核时，我收到不同的结果。如果有人能够发现错误或提供替代解决方案，我将不胜感激！ 这里是我的内核：

__global__ void cuda_dot_kernel(int *n,double *a, double *b, double *dot_res) 
{ 
    __shared__ double cache[threadsPerBlock]; //thread shared memory 
    int global_tid=threadIdx.x + blockIdx.x * blockDim.x; 
    int i=0,cacheIndex=0; 
    double temp = 0; 
    cacheIndex = threadIdx.x; 
    while (global_tid < (*n)) { 
     temp += a[global_tid] * b[global_tid]; 
     global_tid += blockDim.x * gridDim.x; 
    } 
    cache[cacheIndex] = temp; 
    __syncthreads(); 
    for (i=blockDim.x/2; i>0; i>>=1) { 
     if (threadIdx.x < i) { 
      cache[threadIdx.x] += cache[threadIdx.x + i]; 
     } 
     __syncthreads(); 
    } 
    __syncthreads(); 
    if (cacheIndex==0) { 
     *dot_res=cuda_atomicAdd(dot_res,cache[0]); 
    } 
} 

这里是我的设备功能atomicAdd：

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

    return old; 
} 

Answer 1

您正在使用不正确的cuda_atomicAdd功能。内核的这一部分：

if (cacheIndex==0) { 
    *dot_res=cuda_atomicAdd(dot_res,cache[0]); 
} 

是罪魁祸首。在这里，你自动添加到dot_res。然后非原子地设置dot_res与它返回的结果。该函数的返回结果是之前的值，该位置被原子更新，并且它仅提供给“信息”或本地使用调用者。你不会将它分配给你原子更新的东西，这完全违背了使用原子内存访问的目的。代之以这样做：

if (cacheIndex==0) { 
    double result=cuda_atomicAdd(dot_res,cache[0]); 
} 

Answer 2

获取使用ad hoc CUDA代码可能会非常棘手的减少右，所以这里使用推力算法的替代解决方案，它包含在CUDA工具包中：

#include <thrust/inner_product.h> 
#include <thrust/device_ptr.h> 

double do_dot_product(int *n, double *a, double *b) 
{ 
    // wrap raw pointers to device memory with device_ptr 
    thrust::device_ptr<double> d_a(a), d_b(b); 

    // inner_product implements a mathematical dot product 
    return thrust::inner_product(d_a, d_a + n, d_b, 0.0); 
} 

Answer 3

未检查您的代码深度，但这里有一些建议。
如果您仅将GPU用于这样的通用任务，那么我只会建议使用Thrust，因为如果出现复杂的问题，人们不知道要在GPU上高效地编程并行。

1. 启动一个新的并行压缩内核来总结点积。
   由于数据已经存在于设备上，因此启动新内核时性能不会下降。

2. 你的内核看起来并没有在最新GPU上的可能块的最大数量上进行缩放。如果它和你的内核能够计算出数百万个数值的点积，那么由于序列化的原子操作，性能会显着降低。

3. 初学者错误：您的输入数据和共享内存访问范围检查？或者你确定输入数据总是你块大小的倍数？否则你会阅读垃圾。我错误的结果大部分是由于这个错误。

4. 优化您的并行减少。 My Thesis或Optimisations Mark Harris

未经检验的，我只是把它写下来，在记事本：

/* 
* @param inCount_s unsigned long long int Length of both input arrays 
* @param inValues1_g double* First value array 
* @param inValues2_g double* Second value array 
* @param outDots_g double* Output dots of each block, length equals the number of blocks 
*/ 
__global__ void dotProduct(const unsigned long long int inCount_s, 
    const double* inValuesA_g, 
    const double* inValuesB_g, 
    double* outDots_g) 
{ 
    //get unique block index in a possible 3D Grid 
    const unsigned long long int blockId = blockIdx.x //1D 
      + blockIdx.y * gridDim.x //2D 
      + gridDim.x * gridDim.y * blockIdx.z; //3D 


    //block dimension uses only x-coordinate 
    const unsigned long long int tId = blockId * blockDim.x + threadIdx.x; 

    /* 
    * shared value pair products array, where BLOCK_SIZE power of 2 
    * 
    * To improve performance increase its size by multiple of BLOCK_SIZE, so that each threads loads more then 1 element! 
    * (outDots_g length decreases by same factor, and you need to range check and initialize memory) 
    * -> see harris gpu optimisations/parallel reduction slides for more informations. 
    */ 
    __shared__ double dots_s[BLOCK_SIZE]; 


    /* 
    * initialize shared memory array and calculate dot product of two values, 
    * shared memory always needs to be initialized, its never 0 by default, else garbage is read later! 
    */ 
    if(tId < inCount_s) 
     dots_s[threadIdx.x] = inValuesA_g[tId] * inValuesB_g[tId]; 
    else 
     dots_s[threadIdx.x] = 0; 
    __syncthreads(); 

    //do parallel reduction on shared memory array to sum up values 
    reductionAdd(dots_s, dots_s[0]) //see my thesis link 

    //output value 
    if(threadIdx.x == 0) 
     outDots_g[0] = dots_s[0]; 

    //start new parallel reduction kernel to sum up outDots_g! 
} 

编辑：删除不必要的点。