用CUDA C进行点积运算

最近在学习《GPU高性能编程CUDA实战》（机械工业出版社）这本书，但是
在阅读到点积运算这一节时遇到了一些不懂得问题。最近刚刚解决，所以来
与大家分享一下，如果有错误的地方，请大家指点一下。

由于点积运算中核函数是关键也是最难懂的地方，因此在这里我只详细介绍
一下核函数的部分。首先我阐释一下大致的思路。按照书中的示例，进行点
积运算的两个向量长度为33*1024，其中共使用了32个线程块，每个线程块
中使用了256个线程。我们这里也不做改变了。（详情请参考本书第五章内容）

step1:每个线程计算两个元素的积并保存在中间变量temp里（实际计算过程中由于向量长度过长，一次并行计算可能会计算不完，需要迭代几次，因此temp保存的值可能为多个元素积之和），如下图。



相应的代码如下：

__global__ void DotKernel(float *a, float *b, float *c)
{
//声明一个共享内存缓冲区以保存每个线程计算的加和值
//因为编译器将为每个线程块生成一个共享变量的副本
//  所以大小为每个线程块中线程的数量
__shared__ float cache[threadsPerBlock];
int tid = threadIdx.x + blockIdx.x * blockDim.x;
int cacheIndex = threadIdx.x;

float temp = 0;
while (tid < N)
{
temp += a[tid] * b[tid];
tid += blockDim.x * gridDim.x;
}
...
...
}

step2:将temp的值保存到每个线程块共享内存(shared memory)中的数组cache,相应的代码如下：

//设置cache上相应位置上的值
//temp上计算出的值应该为每个线程上几次迭代后计算的加和值
cache[cacheIndex] = temp;

//对线程块中的线程进行同步
//以保证每个线程都完成前面的工作后再执行该语句之后的工作
__syncthreads();

step3:使用归约算法将每个共享内存中的cache归约为一个值，相应的代码如下：

//归约算法将每个线程块上的cache数组归约为一个值cache[0]，最终保存在数组c里
int i = blockDim.x /2;
while (i != 0)
{
if (cacheIndex < i)
cache[cacheIndex] += cache[cacheIndex + i];
__syncthreads();        //确保每个线程已经执行完前面的语句

i /= 2;
}

ps:关于归约算法本书中有详细的介绍。

step4:在每个线程块中选择每个线程索引为0的线程，将归约后的结果传递到全局内存(global memory)中，相关代码如下：

//选择每个线程块中线程索引为0的线程将最终结果传递到全局内存中
if (cacheIndex == 0)
c[blockIdx.x] = cache[0];

到这里核函数的部分已经结束。主函数将调用核函数，并完成剩余的工作。

ps:本博客的内容仅供参考，其中的内容只是为了帮助大家理解《GPU高性能编程CUDA实战》这本书中的内容。