最近在研究并行运算的规约算法，在看《GPU高性能编程CUDA实战》这本书中点积运算时，有些问题想了很久，记录下来；

注点积公式：（dot(A,B)=a1*b1+a2*b2+...+an*bn）

书上例子算点积运算时分为了以下几步：

1.申请共享内存cache;

__shared__ float cache[threadsPerBlock]

这里要注意一个概念：一旦这样声明共享内存，就会创建与线程块的数量相同的数组cache，即每个线程块都会对应一个这样的数组cache且不同线程块中的共享内存是无法交流的;

在这个例子中，共享内存的大小与每个线程块中的线程个数相同；

2.将每个线程块中每个线程计算的元素（a1*b1形式）加起来并放入共享内存；

讯享网int tid = threadIdx.x + blockIdx.x * blockDim.x; int cacheIndex = threadIdx.x; float temp = 0; while (tid < N)//N为数组元素个数 { temp += a[tid] * b[tid]; tid += blockDim.x * gridDim.x; } cache[cacheIndex]=temp; __syncthreads(); 

该例线程块和线程都申请的是一维的，所以tid就代表所有线程块中的所有线程数量；

在while循环这困扰了我很久，不明白为什么需要求和？直到看了博客https://blog.csdn.net/github_/article/details/

原来核函数申请的线程数不一定等于或大于元素个数，也可以小于，这样就需要并行多次，所以这里tid加上了总线程的个数（blockDim.x * gridDim.x）；

3.将每个线程块中的每个线程的元素加起来；

int i = blockDim.x /2;//单个线程块中线程数的一半 while (i != 0) { if (cacheIndex < i) cache[cacheIndex] += cache[cacheIndex + i]; __syncthreads(); i /= 2; } 

此处就用到了并行规约的思想，规约书上的定义是这样的：对一个输入数组执行某种计算，然后产生一个更小的结果数组，这个过程就叫规约；；

https://www.cnblogs.com/5long/p/algorithms-on-cuda-reduction.html这个博客规约讲得挺好；

讯享网

本例此处的思想就是最右侧的图示，通过并行的方法求两两元素之和；

例如假设本例 blockDim.x=8，那么i=4;

执行第一次并行就为[0,4],[1,5],[2,6],[3,7]，超过或等于i的索引就不计算了;

然后i=2;计算[0,2],[1,3]

...

直到i==0;

最后的值就保存在了cache[0]中;然后把cache[0]给全局变量c后，我们就得到了每个线程块所有线程的元素和；

讯享网 if (cacheIndex == 0) c[blockIdx.x] = cache[0];

这里cacheIndex可以为不超过线程索引的任何数，因为所有的线程都执行一样的操作，因此一个就行；

4.将每个线程块中的线程元素和加起来得到最终结果

这一步已经不在核函数了，例子是将数组拷贝出来在cpu中计算出来的，这里很简单，不再赘述了。

书中完整代码：

#include<iostream> #include"cuda_runtime.h" #define min(a,b) (a<b ? a:b) #define sum_squares(x) (x*(x+1)*(2*x+1)/6) const int N=33*1024; const int threadsPerBlock=256; const int blockPerGrid=min(32,(N+threadsPerBlock-1)/threadsPerBlock); using namespace std; __global__ void dot(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; } cache[cacheIndex]=temp; __syncthreads(); //对于归约运算来说，以下代码要求threadPerBLock必须为2的指数 int i=blockDim.x/2; while(i != 0) { if(cacheIndex<i) cache[cacheIndex] += cache[cacheIndex+i]; __syncthreads(); i /=2; } if(cacheIndex==0) { c[blockIdx.x]=cache[0]; } } int main() { float *a,*b,c,*partial_c; float *dev_a,*dev_b,*dev_partial_c; a=(float*)malloc(N*sizeof(float)); b=(float*)malloc(N*sizeof(float)); partial_c=(float*)malloc(blockPerGrid*sizeof(float)); cudaMalloc((void)&dev_a,N*sizeof(float)); cudaMalloc((void)&dev_b,N*sizeof(float)); cudaMalloc((void)&dev_partial_c,blockPerGrid*sizeof(float)); for(int i=0;i<N;i++) { a[i]=i; b[i]=i*2; } cudaMemcpy(dev_a,a,N*sizeof(float),cudaMemcpyHostToDevice); cudaMemcpy(dev_b,b,N*sizeof(float),cudaMemcpyHostToDevice); dot<<<blockPerGrid,threadsPerBlock>>>(dev_a,dev_b,dev_partial_c); cudaMemcpy(partial_c,dev_partial_c,blockPerGrid*sizeof(float),cudaMemcpyDeviceToHost); c=0; for(int i=0;i<blockPerGrid;i++) { c+=partial_c[i]; } printf("c = %f\n",c); cudaFree(dev_a); cudaFree(dev_b); cudaFree(dev_partial_c); free(a); free(b); free(partial_c); }