[CUDA学习]4.线程协作

线程协作

目标

• 了解CUDA C中的线程
• 了解不同线程之间的通信机制
• 了解并执行线程的同步机制

并行线程块的分解

前面我们了解了如何在GPU上启动并行代码。具体的方法是，告诉CUDA运行时启动核函数的多个并行副本。这里将这些并行副本称为线程块(Block)。
CUDA运行时将这些线程块分解为多个线程，当需要启动多个并行线程块时，只需要将尖括号中的第一个参数由1改为想要启动的线程块的数量。例如在介绍矢量加法示例的时候，曾经就是为矢量中的每个元素(共有N个元素)都启动一个线程块。调用如下：

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

在尖括号中，第二个参数表示CUDA运行时在每个线程块中创建的线程数量。上面这行代码表示每个线程块中只启动一个线程。总共启动线程的数量可以按下面公式计算：
N个线程块*1个线程/线程=N个并行线程
事实上我们也可以启动N/2个线程块，每个线程块包含两个线程，或者N/4个线程块，每个线程块包含4个线程。

矢量求和：重新回顾

这里将实现一个和前面相同的任务，即将两个输入矢量相加并将得到的结果保存在第三个输出矢量中。这里将使用线程而不是线程块来实现这一工作。这两种方法现阶段没有任何区别，但是，线程块中的并行线程能够完成并行线程块无法完成的工作。

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

在将实现方式由并行线程块改为并行线程时，首先需要修改两个地方。第一个地方是，之前核函数的调用中是启动N个线程块，并且每个线程块对应一个线程，具体如下所示：

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

现在我们将启动N个线程，所有的线程都在一个线程块中：

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

第二个要修改的地方就是对数据的索引的方法。之前核函数中是通过线程块索引对输入数据和输出数据进行索引。

int tid = blockIdx.x;

这里由于只有一个线程块，所以将采用线程索引对数据进行索引。

int tid = threadIdx.x;

所以采用并行线程实现GPU上的矢量求和的例子：

#include <stdio.h>#define N 10__global__ void add(int *a,int *b,int *c){    int tid = threadIdx.x;    if(tid < N){        c[tid] = a[tid] + b[tid];     }}int main(void){    int a[N],b[N],c[N];    int *dev_a,*dev_b,*dev_c;    //在GPU上分配内存    cudaMalloc((void**)&dev_a,N*sizeof(int));    cudaMalloc((void**)&dev_b,N*sizeof(int));    cudaMalloc((void**)&dev_c,N*sizeof(int));    //在CPU上为数组‘a’和数组‘b’赋值    for(int i=0;i<N;i++){        a[i] = i;        b[i] = i*i;    }    //将数组‘a’和数组‘b’复制到GPU上    cudaMemcpy(dev_a,a,N*sizeof(int),cudaMemcpyHostToDevice);    cudaMemcpy(dev_b,b,N*sizeof(int),cudaMemcpyHostToDevice);    add<<<1,N>>>(dev_a,dev_b,dev_c);    //将数组‘C’从GPU复制到CPU中    cudaMemcpy(c,dev_c,N*sizeof(int),cudaMemcpyDeviceToHost);    //显示结果    for(int i=0;i<N;i++){        printf("%d+%d=%d\n",a[i],b[i],c[i]);    }    //释放内存    cudaFree(&dev_a);    cudaFree(&dev_b);    cudaFree(&dev_c);}

在GPU上求更长矢量的和

在前面我们了解到硬件将线程块的数量限制为不超过65535个。同样对于启动核函数时每个线程块中的线程数量，硬件也做了限制。具体的来说，最大线程数量不能超过设备属性中的maxThreadsPerBlock域的值，所以对于一些比较大的矩阵就需要将线程与线程块结合起来实现计算。这种方法主要需要修改两个地方：
* 核函数的索引计算方法int tid = threadIdx.x + blockIdx.x*blockDim.x这里面使用了一个新的变量blockDim。对于所有的线程来说这个变量是一个常数，保存的是线程块中每一维的线程的数量。由于使用的是使用一维线程线程块，因此只需要用到blockDim.x。回顾第四章的内容，在gridDim中保存了一个类似的值，即在线程格每一维的线程块数量。此外gridDim是二维的，而blockDim实际上是三维的。也就是说CUDA运行时允许启动一个二维的线程格，并且线程中每个线程块都是一个三维的线程数组。
* 核函数的调用方式。这种方法需要修改核函数调用本身。虽然任然需要N个并行线程，但是我们希望在多个线程块中启动这些线程，这样就不会出现每个线程块最多不会超过512个线程的最大数量的限制。这里将每个线程块中包含的线程的数量设置为128。然后，可以启动N/128个线程块，这样总共就启动了N个线程同时运行。这里的问题在于，N/128是一个整除，所以这里我们希望能够向上取整数。具体核函数调用方法如下：add<<<(N+127)/128,128>>>(dev_a,dev_b,dev_c)。这里会发现将启动过多的线程。然而，在核函数中已经解决了这一问题：

if(tid<N){    c[tid] = a[tid]+b[tid];}

因此，当索引越过数组边界时，例如当启动的并行线程数量不是128的整数倍的时候就会出现这种情况，那么核函数将自动停止计算。更重要的是，核函数不会对越过数组边界的内存进行存储或者读写操作。

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

当第一次介绍在GPU上启动并行线程块的时候，我们并没有考虑所有可能的情况，除了在线程数量商存在限制之外，在线程块的数量上同样存在着一个硬件的限制，线程格的每一个维度的大小不能超过65535。
因此，这就对矢量加法的实现带来一个问题。如果启动N/128个线程相加，那么当矢量长度超过65535 X 128的时候核函数就会调用失败。这里需要对核函数做下面的修改：

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

这里可以看出当使用的线程数量未达到要计算的量，该线程不会在进行下一个计算。在GPU系统中我们将并行线程的数量看作是处理器的数量。尽管GPU处理单元的数量可能小于这个值，但是我们认为每个线程在逻辑单元上都可以并行执行，并且硬件可以调度这个线程以便实际执行。通过并行化过程与硬件的实际执行过程解耦开来。
在理解了上述实现方式所包含的核心思想后，我们只需要知道如何计算每个并行线程的初始索引，以确定递增量的值。我们希望每个并行线程从不同的索引开始，因此就需要对线程索引和线程块索引进行线性化。int tid = threadIdx.x + blockIdx.x*blockDim.x。在完成上面的计算后就需要对索引的增量进行计算，这个数值等于每个线程块中线程的数量乘以线程格中线程块的数量。
在核函数调用的时候为了保证不启动过多的线程块，我们将线程固定为某个较小的值。完整的代码如下：（这里计算的矢量长度取决与GPU上的内存容量）

#include <stdio.h>#define N (33*1024)__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;    }}int main(void){    int a[N],b[N],c[N];    int *dev_a,*dev_b,*dev_c;    //在GPU上分配内存    cudaMalloc((void**)&dev_a,N*sizeof(int));    cudaMalloc((void**)&dev_b,N*sizeof(int));    cudaMalloc((void**)&dev_c,N*sizeof(int));    //在CPU上为数组'a'和数组'b'赋值    for(int i=0;i<N;i++){        a[i] = i;        b[i] = i*i;    }    //将数组‘a’和数组‘b’复制到GPU内存中    cudaMemcpy(dev_a,a,N*sizeof(int),cudaMemcpyHostToDevice);    cudaMemcpy(dev_b,b,N*sizeof(int),cudaMemcpyHostToDevice);    add<<<128,128>>>(dev_a,dev_b,dev_c);    //将数组‘C’从GPU复制到CPU中    cudaMemcpy(c,dev_c,N*sizeof(int),cudaMemcpyDeviceToHost);    //验证GPU计算结果    for(int i=0;i<N;i++){        if(a[i]+b[i] != c[i]){            printf("error:%d+%d=%d\n",a[i],b[i],c[i]);        }    }    cudaFree(dev_a);    cudaFree(dev_b);    cudaFree(dev_c);    return 0;}

在GPU上使用线程实现纹波效果

#include "cuda.h"#include "../common/book.h"#include "../common/cpu_anim.h"#define DIM 1024#define PI 3.1415926535897932/* 核函数需要知道当前动画时间以便于生成正确的帧。在CPUAnimBitmap的代码中 * 将当前时间ticks传递给generate_frame() * /__global__ void kernel(unsigned char *ptr,int ticks){    int x = threadIdx.x + blockIdx.x*blockDim.x;    int y = threadIdx.y + blockIdx.y*blockDim.y;    int offset = x + y*blockDim.x*gridDim.x;    float fx = x - DIM/2;    float fy = y - DIM/2;    float d = sqrtf(fx*fx + fy*fy);    unsigned char gray = (unsigned char)(128.0f + 127.0f*                                                        cos(d/10.0f - ticks/7.0f)/                                                        (d/10.0f + 1.0f));    ptr[offset*4 + 0] = gray;    ptr[offset*4 + 1] = gray;    ptr[offset*4 + 2] = gray;    ptr[offset*4 + 3] = 255;}struct DataBlock{    unsigned char *dev_bitmap;    CPUAnimBitmap *bitmap;};void generate_frame(DataBlock *d,int ticks){    dim3 blocks(DIM/16,DIM/16);    dim3 threads(16,16);    kernel<<<blocks,threads>>>(d->dev_bitmap,ticks);    cudaMemcpy(d->bitmap->get_ptr(),               d->dev_bitmap,               d->bitmap->image_size(),               cudaMemcpyDeviceToHost);}//释放在GPU上分配的内存void cleanup(DataBlock *d){    cudaFree(d->dev_bitmap);}int main(void){    DataBlock data;    CPUAnimBitmap bitmap(DIM,DIM,&data);    data.bitmap = &bitmap;    cudaMalloc((void**)&data.dev_bitmap,bitmap.image_size());    /*     * 这里将一个指向generate_frame()的函数指针传递给anim_and_exit()。每当要生成一副新的动画     * 时都要调用函数generate_frame();     */    bitmap.anim_and_exit((void(*)(void*,int))generate_frame,//函数的指针调用                            (void(*)(void*))cleanup);    return 0;}

共享内存和同步

到目前为止将线程块分解为线程的目的只是为了解决线程块数量的硬件限制。这是一个很勉强的动机，因为这个完全可以由CUDA运行时在幕后完成。显然还有其他原因需要将线程分解为多个线程。
CUDA C支持共享内存。这块内存引出了对C语言的另一个扩展，这个扩展类似于__device__和__global__。在编写代码的时候，你可以将CUDA C的关键字__share__添加到变量的声明中，使得这个变量存储在共享内存中。那么这样做的目的是什么了？
CUDA C编译器对共享内存中的变量和普通变量将分别采取不同的处理方式。对于在GPU上自动启动的多个线程块，CUDA C编译器都将创建该变量的一个副本。线程块中的每个线程都将共享这块内存，但是线程中无法看到也不能修改其他线程块变量副本。这就实现了一种非常好的方式，使得一个线程块中的多个线程能够在计算上进行通信和协作。而且共享内存的缓存区驻留在物理GPU上，而不是驻留在GPU之外的系统内存中。因此，在访问共享内存时的延迟要远远低于访问普通缓冲区的延迟，使得共享内存像在每个线程块的高速缓存或者中间结果暂存器上那样高效。
线程之间进行通讯，那么还需要一种机制来实现县城之间的同步。例如：如果线程A将一个值写入到共享内存，并且我们希望线程B对这个值进行一些操作，那么只有当线程A将写入操作完成后，线程B才能开始执行它的操作。如果没有同步，那么将发生竞态条件，在这种情况下，代码执行结果的正确性将取决硬件的不确定性。

点积运算

矢量的点积运算计算公式如下：
(x1,x2,x3,x4)·(y1,y2,y3,y4)=x1y1+x2y2+x3y3+x4y4
这里可以像矢量加法那样，每个线程将两个相应的元素相乘，然后移动到下两个元素。由于最终的结果是所有乘积的总和，因此每个线程还要保存它所计算的成积的加和。与矢量加法类似，线程每次对索引的增加值为线程的数量，这样能确保不会遗漏任何元素。而且也不会将任何元素相乘两次。具体代码如下：

#include <stdio.h>#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){    /*     *  代码在这里声明了一个内存共享缓冲区，名字为cache，这个缓存区保存每个     *  线程计算的加和值，这里将数组大小声明为threadsPerBlock，这样线程块中的     *  每个线程都能被分配足够的内存，即线程块的数量乘以threadsPerBlock。但对     *  于共享变量来说，由于编译器将为每个线程块生成一个共享变量的副本，因此     *  只需要根据线程块中线程的数量来分配内存。     *       */    __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    }    //设置catch中相应位置的值。这里设置了共享内存缓冲区相应位置上的值，    //以便随后能将整个数组相加起来，并无需担心某个特定位置商是否包含有效的数据    catche[cacheIndex] = temp;    /*     * 当算法执行到这一步时，我们需要将cache中所有的值相加起来。在执行这个运算时，     * 需要通过一个线程来读取保存在cache中的值，然而前面提到过，这可能是一种危险的     * 操作。我们需要某种方法来确保所有对共享数组cache[]的写入操作在读取的cache之前     * 完成。具体为使用__syncthreads();对线程同步。这个函数将确保线程块中的每个线程都     * 执行完__syncthreads()前面的语句，才会执行后面的第一条语句，这时候，线程块中的     * 其他线程肯定都执行完了__syncthreads()。     */     //对线程块同步     __syncthreads();     //通过归约运算，计算每个线程块内的和，这里要求threadPerBlock必须是2的指数     int i = blockDim.x/2;     while(i != 0){        if(cacheIndex<i){            catche[catcheIndex] += catche[catcheIndex + i];            __syncthreads();            i /= 2；        }        if(cacheIndex == 0)            c[blockIdx.x] = cache[0];     }}int main(void){    float *a, *b, c, pratial_c;    float *dev_a,*dev_b,*dev_pratial_c;    //在CPU上为数组分配内存    a = (float*)malloc(N*sizeof(float));    b = (float*)malloc(N*sizeof(float));    partial_c = (float*)malloc(blocksPerGrid*sizeof(float));    //在GPU上分配内存    cudaMalloc((void**)&dev_a,N*sizeof(float));    cudaMalloc((void**)&dev_b,N*sizeof(float));    cudaMalloc((void**)&dev_pratial_c,blocksPerGrid*sizeof(float));    //填充主机内存    for(int i=0;i< N;i++){        a[i] = i;        b[i] = i*2;    }    //将数组‘a’和数组‘b’复制到GPU上    cudaMemcpy(dev_a,a,N*sizeof(float),cudaMemcpyHostToDevice);    cudaMemcpy(dev_b,b,N*sizeof(float),cudamemcpyHostToDevice);    dot<<<blocksPerGrid,threadsPerBlock>>>(dev_a,dev_b,dev_partial_c);    //将数组‘C’从GPU复制到CPU    cudaMemcpy(partial_c,dev_partial_c,blocksPerGrid*sizeof(float),cudaMemcpyDeviceToHost);    //在CPU上完成最后的求和运算    c = 0;    for(int i = 0;i<blocksPerGrid;i++){        c += partial_c[i];    }    #define sum_squares(x) (x*(x+1)*(2*x+1)/6)    printf("Does GPU Value %.6g = %.6g?\n",c,2*sum_squares((float)(N-1)));    //释放GPU上的内存    cudaFree(dev_a);    cudaFree(dev_b);    cudaFree(dev_partial_c);    //释放CPU上的内存    free(a);    free(b);    free(partial_c);}

（错误的）点积运算优化

在前面介绍了点积运算示例中的第二个__syncthreads()。现在我们将在这里进一步分析这个函数的调用，并尝试对这个函数做出改进。之所以需要第二个__syncthreads()，是因为在循环迭代中更新了共享内存变量cache()，并且在循环的下一次迭代开始之前，需要确保当前迭代中所有线程的更新操作都已经完成。

int i = blockDim.x/2;while(i != 0){    if(cacheIndex < i){        cache[cacheIndex] += cache[cacheIndex + i];    }    __syncthreads();    i /= 2;}

在上面的程序中，我们发现只有当cacheIndex小于i的时候才需要更新共享内存的缓存区cache[]。由于cacheIndex实际上就是等于threadIdx.x，因此这意味着只有一部分线程会更新共享内存。由于调用__syncthreads的目的只是为了确保这些更新在下一次迭代之前已经完成，因此如果将代码修改为了等待那些需要写入共享内存的线程，是不是就能获得性能提升了。

int i = blockDim.x/2;while(i != 0){    if(cacheIndex < i){        cache[cacheIndex] += cache[cacheIndex + i];         __syncthreads();    }    i /= 2;}

虽然这种优化代码的初衷不错，但是却不起作用。实际上这种情况比优化之前的情况还要糟糕。再这样对核函数修改之后，GPU将停止响应，因此不得不强行终止程序。在GPU中的线程块中的每个线程依次通过代码，每次一行。每个线程都要执行相同的命令，但是对不同的数据进行运算。然而，当每个线程执行的指令位于一个条件语句中，如if()语句中的时候，那么将会出现什么情况。显然这时候不是每个线程都会执行这个指令。当某一些线程需要执行一条指令，而其他线程不需要执行时，这种情况就称为线程的发散。在__syncthreads()情况下。应为CUDA架构将确保，除非线程块中的每个线程都执行__syncthreads()，否则没有任何线程能执行__syncthreads()之后的指令。遗憾的是。如果__syncthreads()处于发散分支中，那么一些线程将永远无法执行__syncthreads()。因此由于要确保每个线程执行完__syncthreads()后才能执行后面的语句，因此硬件将使这些线程保持等待。所以。如果在点积运算示例中将__syncthreads()放到if语句中，那么任何cacheIndex大于或等于i的线程将永远不能执行__syncthreads()。这将使处理器挂起，应为GPU在等待某个永远不可能发生的事件。

基于共享内存的位图

#include "cuda.h"#include "../common/book.h"#include "../common/cpu_bitmap.h"#define DIM 1024#define PI 3.1415926535897932f__global__ void kernel( unsigned char *ptr ) {    // map from threadIdx/BlockIdx to pixel position    int x = threadIdx.x + blockIdx.x * blockDim.x;    int y = threadIdx.y + blockIdx.y * blockDim.y;    int offset = x + y * blockDim.x * gridDim.x;    __shared__ float    shared[16][16];    // now calculate the value at that position    const float period = 128.0f;    shared[threadIdx.x][threadIdx.y] =            255 * (sinf(x*2.0f*PI/ period) + 1.0f) *                  (sinf(y*2.0f*PI/ period) + 1.0f) / 4.0f;    //由于最后是通过逆序复制的不加线程同步会导致未被计算先给了一个随机的值    __syncthreads();    ptr[offset*4 + 0] = 0;    ptr[offset*4 + 1] = shared[15-threadIdx.x][15-threadIdx.y];    ptr[offset*4 + 2] = 0;    ptr[offset*4 + 3] = 255;}// globals needed by the update routinestruct DataBlock {    unsigned char   *dev_bitmap;};int main( void ) {    DataBlock   data;    CPUBitmap bitmap( DIM, DIM, &data );    unsigned char    *dev_bitmap;    HANDLE_ERROR( cudaMalloc( (void**)&dev_bitmap,                              bitmap.image_size() ) );    data.dev_bitmap = dev_bitmap;    dim3    grids(DIM/16,DIM/16);    dim3    threads(16,16);    kernel<<<grids,threads>>>( dev_bitmap );    HANDLE_ERROR( cudaMemcpy( bitmap.get_ptr(), dev_bitmap,                              bitmap.image_size(),                              cudaMemcpyDeviceToHost ) );    HANDLE_ERROR( cudaFree( dev_bitmap ) );    bitmap.display_and_exit();}

本章小结

这里了解了如何将线程块进一步分解为更小的并行执行单元。这种并行执行单元也称为线程。我们回顾前面矢量相加的示例，并介绍如何实现任意长度矢量的加法。我们还给出了一个规约运算的示例，以及如何通过 共享内存和同步来实现这个运算。这个示例说明了如何对GPU和CPU进行协作完成运算。最后，我们还给出了当忽略同步的时候给应用程序造成的问题。
到这里你已经学习了CUDA C的大部分基础概念，它与标准C在某些方面是相同的，但在大多数方面是不同的。此时，你可以思考你曾经遇到过的问题，以及那些问题可以通过CUDA C的并行实现。后面将讲解CUDA C的更多功能以及一些高级的原子函数。