最近在研究并行运算的规约算法，在看《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_30605157/article/details/48196997
原来核函数申请的线程数不一定等于或大于元素个数，也可以小于，这样就需要并行多次，所以这里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);
}
