CUDA基础的基础教程二、线程协作

1.回顾：

上一篇博客我们介绍了CUDA调用核函数的形式 <<<N,M>>>> 。
N 代表线程块的数量，M 代表线程块中有多少个并行的线程，那么就是创建了N×M个并行线程，当然 M=1 时就是 N 个，我们利用 CUDA 这个性质，重新回顾下上一篇说到的矢量求和，我们之前用的是线程块，我们用线程来写的话，只要变成<<<1,N>>> 就行了，但是区别在哪里呢？
由于每个 GPU 能支持的最大进程块的数量不能超过 65535，每个块内的线程数不能超过 512，注意每个硬件的规格可能是不一样的。于是我们只能将这两个东西结合一下了。其他的都不变，只需要改变一下定位的方式就可以了

2.任意长度的矢量求和：

首先我们要对线程块再度包装，介绍一下线程格。线程格相当于是一个二维的平面，每个平面上的格子都是一个线程块，这样就方便找到了线程块。

线程格需要预先定义 dim3 gird(DIM,DIM) ， dim3 代表一个三维数组（封装的）但是一般来说不声明第3维的话就是1，这个语句声明了一个 DIM×DIM 大小的线程格，调用的时候，要用dim3这个变量代替第一个参数。

我们想要知道每个线程块的位置，就可以使用预先定义的 gridIdx.x 代表这个格子的“x坐标”也就是第一维的位置。进程格的大小就是 gridDim.x,gridDim.y

同理，我们想知道某个线程的位置，就需要用到 threadIdx.x 这个变量，我们来列个表格

位置参数	线程格	线程块	线程
大小	gridDim	blockDim	…
位置	…	blockIdx	threadIdx

我们知道了这些预先定义的变量，就知道怎么找到正确的位置了

3：改装函数

我们先上改装过的函数：

__global__ void add(int *a,int *b,int *c)
{
    int tid = threadIdx.x + blockIdx.x * blockDim.x;
    while(tid < N) 
    {
        c[tid] = a[tid] + b[tid];
        tid += blockDim.x * gridDim.x;
    }
}

首先我们这里进行了重新的取值，新的 tid 就是这个线程在相关的线程块的位置，然后就是确定增量，我们给每个进程多分配一点任务，每次给tid一个增量，这个增量是我们现在正在并行的线程数,这个地方不太好理解，画画图就ok了.
其他和上一个例子一样就ok，但是调用的时候，我们要使用<<<128,128>>>，或者根据具体的大小进行设置

4. 一个新的例子：点积运算

两个向量的点积是:
$$(x1,x2,x3,x4)\ast (y1,y2,y3,y4)=x1y1+x2y2+x3y3+x4y4$$
我们使用CUDA进行优化，首先还是写个核函数，将问题分解成对应位置的相乘，然后加在一起。我们首先完成核函数的划分：
我们首先需要每个进程计算两个对应位置的乘积，然后顺次移动到下一个位置，移动的增量是线程的数量

	__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;

这里先解释一下 __shared__ 的作用，这里这个修饰符修饰了一个数组，代表了正式个共享的内存空间，相当于在这个数组是一个静态的变量（CPU下）。他的大小是每个线程快的线程数，也就是我们给一个线程块内的线程分配了一个全局的变量，每个线程将自己的答案存到正确的位置上，然后我们就要将这个线程块的答案累加到一起。进行这一步之前，我们要确保每个线程已经结束了第一步的乘法运算，因为我们要对共享内存进行读取，如果没有计算完成，就会出现问题。

我们调用 __syncthreads() 运行之前的代码之后在这里等待其他的线程也运行到此。
我们现在需要对共享数组进行一个规约 (Reduction) 操作，也就是求和的fancy说法，我们分析一下，会发现这个操作也是可以并行的，那我们就利用已经有的线程继续计算：

     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];
 

这里每个线程每次只选择两个数操作，然后我们等待同步，这时只有前一半线程参与了工作，并且重写了前一半的内存，然后我们参与线程的数字减半，然后只对刚才一半被重写的内存接着进行重写……依次类推，最终答案会被保存在 cache[0] 中。
这里还有一个优化的地方，我们只让第0个线程进行了写入操作，实际上，我们每个线程都可以进行写入操作，但是我们何必产生那么多信号量呢？

下面是核函数完整的代码，之后在CPU上进行一步求和就行了，因为每个线程块分配到的任务是不一样的，但是最麻烦的操作已经完成，剩下的放到CPU上就可以了。

__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;
     }
     
     // set the cache values
     cache[cacheIndex] = temp;
     
     // synchronize threads in this block
     __syncthreads();
 
     // for reductions, threadsPerBlock must be a power of 2
     // because of the following code
     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];
 }