GPU 高性能编程 CUDA : 线程协作

并行线程块的分解：

在矢量加法中，为矢量中的每一个元素都启动一个线程块

add<<<N,1>>>(dev_a,dev_b,dev_c);

尖括号中的第一个参数创建的线程块的数量，第二个参数表示每个线程块中创建的线程数量，所以上述启动的线程数量为 N （     N*1）

使用线程实现 GPU 上的矢量求和：

需要修改两个地方，第一个地方是

add<<<N,1>>>(dev_a,dev_b,dev_c);
add<<<1,N>>>(dev_a,dev_b,dev_c);

第二个地方是

int tid=blockIdx.x;
int tid=threadIdx.x;

硬件线程块的数量限制为 65535，每个线程块中的线程数量不超过设备属性中的 maxTreadsPerBlock 域的值

在 GPU 上对更长的矢量求和：

和之前一样，改动两个地方：核函数中索引的计算方法和核函数的调用方式

int tid = threadIdx.x + blockIdx.x*blockDim.x;
add<<<(N+127)/128,128>>>(dev_a,dev_b,dev_c);

类似于Julia 中将二维图像索引线性化

示例中假定每一个线程块有 128 个线程数，可以启动 N/128 个线程块进行计算，但是这是一个整除除法，所以加 127 保证可以启动足够多的线程块，在核函数中访问数组合输出数组之前一定要检查索引是否在 0 到 N 之间，保证索引不超过数组的边界

在 GPU 上对任意长度的矢量求和：

如果当矢量的长度大于 65535*129=8388480 时，核函数会调用失败

首先对核函数进行一些修改

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

计算每个并行线程的初始索引和确认递增的量值，先对线程索引和线程块索引进行线性化，之后每个线程对当前索引上的任务计算完成后，对索引进行递增，递增的步长为线程格中正在运行的线程数量，等于每个线程块的线程数量乘以线程格中线程块的数量，即 blockDim.x * gridDim.x

之后为了确保不启动过多的线程块，将线程块的数量定为一个任意固定的值，只要在限制范围内，例如128

add<<<128,128>>>(dev_a,dev_b,dev_c);

共享内存和协作：

CUDA C 编译器对共享内存中的变量和普通变量采取不同的处理方式，对于在 GPU 上启动的每一个线程块， CUDA C 编译器都将创建该变量的一个副本，线程块的每一个线程都共享这块内存，但线程却无法看到也不能修改其他线程块的变量副本。

int tid = threadIdx.x + blockDim.x * blockIdx.x;
int cacheIndex=threadIdx.x;

通过线程索引和线程块索引计算出输入数组中的一个全局偏移，共享内存偏移就等于线程索引，线程块索引和这个偏移无关，因为每个线程块都有该共享内存的私有副本

//对线程块中的线程进行同步
__syncthreads();

这个语句保证每个线程都执行完 __syncthreads() 前面的语句后，才会执行下一句

点击运算：

#define imin(a,b) (a < b ? a: b);
const int N = 33*1024;
const int threadsPerBlock = 256;
const int blocksPerGrid = imin(32,(N+threadsPerBlock-1)/threadsPerBlock);

__global__ void dot(float *a,float *b,float *c){
    __shared__ float cache[threadsPerBlock];
    int tid = threadsIdx.x + BlockDim.x * BlockIdx.x;
    int cacheIndex = threadsIdx.x;

    float temp = 0;

    while(tid < N){
        temp += a[tid]*b[tid];
        tid += BlockDim.x * gridDim.x;
    }
    //设置 cache 中相应位置的值
    cache[cacheIndex]= tmp ;
    //对线程块中的线程进行同步
    _syncthreads();
    //对 cache 内数据进行规约计算
    int i = threadsPerBlock/2;

    while(i !=0){
        cache[threadIdx.x]+=cache[threadIdx.x+i];

        _syncthreads();
        i/=2;
    }

    if(cacheIndex == 0){
        c[blockIdx.x]=cache[0]}
    }}
int main(void){
    float *a,*b,c,*partial_c;
    float *dev_a,*dev_b,*partial_c;

    //CPU上分配内存
    a=(float*)malloc(N*sizeof(float));
    b=(float*)malloc(N*sizeof(flota));
    partial_c=(float*)malloc(blocksPerGrid*sizeof(float));

    //GPU 上分配内存
    cudaMalloc((void**)&dev_a,N*sizeof(float));
    cudaMalloc((void**)&dev_b,N*sizeof(float));
    cudaMalloc((void**)&partial_c,BlocksPerGrid*sizeof(float));
    //填充主机内存
    for（int i=0；i<N;i++）{
        a[i]=i;
        b[i]=i*i;}
    //复制数据到社备
    cudaMemcpy(dev_a,a,N*sizeof(float));
    cudaMemcpy(dev_b,b,N*sizeof(float));

    dot<<<blocksPerGrid,threadPerBlock>>>(dev_a,dev_b,dev_partial_c);
    //将数组 c 复制到 CPU
    cudaMemcpy(partial_c,dev_partial_c,blockPerGrid*sizeof(float),cudaMemcpyDeviceToHost);
    //在 CPU 上完成最后的求和运算
    c = 0;
    for(int i = 0;i<blockPerGrid;i++){
        c+=parital_c[i];
    }
    //释放 GPU 上的内存
    cudaFree(dev_a);
    cudaFree(dev_b);
    cudaFree(dev_partial_c);

    //释放 CPU 上的内存
    free(a);
    free(b);
    free(parital_c);
}

使用共享内存和 __synscthreads 来确保数据在继续进行之前就已经就绪