11.理解线程束执行的本质（一）

  前面已经大概的介绍了CUDA执行模型的大概过程，硬件的大概结构，例如SM的大概结构，而对于硬件来说，CUDA执行的实质是线程束的执行。从软件的角度看，cuda执行中似乎所有的线程都是并行地运行的。在逻辑上这是正确的，但从硬件的角度来看，实际上硬件资源是有限的，不是所有线程在物理上都可以同时并行地执行。现在从硬件的角度来介绍线程束执行，了解执行的硬件过程是有必要的。

1.线程束和线程块

  线程束是SM中基本的执行单元。当一个网格被启动（网格被启动，等价于一个内核被启动，每个内核对应于自己的网格），网格中包含线程块，线程块被分配到某一个SM上以后，将分为多个线程束，每个线程束一般是32个线程，在一个线程束中，所有的线程按照单指令多线程（SIMT）方式执行；也就是说，所有线程都执行相同的指令，每个线程在私有数据上进行操作。下图展示了线程块的逻辑视图和硬件视图之间的关系。

• 逻辑视图：即我们软件的逻辑视角，我们是通过启动核函数来调用cuda执行并行程序的，而我们启动核函数，只指定了线程格和线程块的数量。
• 硬件视图：即硬件根据软件的指令执行的视角，硬件是已经设定好的，不需要我们过多操心的，但理解它的执行过程有助于我们更好的控制cuda.
  
    在计算机里，内存总是一维线性存在的，所以程序执行起来也是一维的访问线程块中的线程，但是我们在写程序的时候却可以以二维三维的方式进行，原因是方便我们写程序。在一个线程块中，每个线程都有一个唯一的ID。例如，一个有128个线程的一维线程块被组织到4个线程里，如下所示。
  
    对于一维的线程块，唯一的线程ID被存储在CUDA的内置变量threadIdx.x中当。当线程块使用三维编号时，x位于最内层，y位于中层，z位于最外层，所以对于一个给定的二维线程块，在一个块中每个线程的独特标识符都可以用内置变量threadIdx和blockDim来计算：

threadIdx.y * blockDim.x + threadIdx.x

  对于一个三维线程块，计算如下：

threadIdx.z * blockDim.y * blockDim.x + threadIdx.y * blockDim.x + threadIdx.x

  一个线程块的线程束的数量可以根据下式确定：
$$\text{一个线程块中线程束的数量}=\text{ceil}\left(\begin{array}{c}\frac{\text{一个线程块中线程的数量}}{\text{线程束大小}}\end{array}\right)$$

ceil函数：是向正无穷取整的函数，即向上取整。

  因此，硬件总是给一个线程块分配一定数量的线程束。线程束不会在不同的线程块之间分离。如果线程块的大小不是线程束大小的偶数倍，那么在最后的线程束里有些线程就不会活跃。例如下图是一个在x轴中有40个线程、在y轴中有2个线程的二维线程块。从应用程序的角度来看，在一个二维网格中共有80个线程。硬件为这个线程块配置了3个线程束，使总共96个硬件线程去支持80个软件线程。注意，最后半个线程束是不活跃的。即使这些线程未被使用，它们仍然消耗SM的资源，如寄存器等。

1.1 线程块：逻辑角度与硬件角度

• 从逻辑角度来看：线程块是线程的集合，它们可以被组织为一维、二维或三维布局。
• 从硬件角度来看：线程块是一维线程束的集合。在线程块中线程被组织成一维布局，每32个连续线程组成一个线程束。

2.线程束分化

  控制流是高级编程语言的基本构造中的一种。GPU支持传统的、C风格的、显式的控制流结构，例如，if…then…else、for和while。

  CPU拥有复杂的硬件以执行分支预测，也就是在每个条件检查中预测应用程序的控制流会使用哪个分支。如果预测正确，CPU中的分支只需付出很小的性能代价。如果预测不正确，CPU可能会停止运行很多个周期，因为指令流水线被清空了。

  GPU是相对简单的设备，它没有复杂的分支预测机制。一个线程束中的所有线程在同一周期中必须执行相同的指令，如果一个线程执行一条指令，那么线程束中的所有线程都必须执行该指令。如果在同一线程束中的线程使用不同的路径通过同一个应用程序，这可能会产生问题。例如，思考下面的语句：

if (con)
{
    //do something
}
else
{
    //do something
}

  假设这段代码是核函数的一部分，那么当一个线程束的32个线程执行这段代码的时候，如果其中16个执行if中的代码段，而另外16个执行else中的代码块，同一个线程束中的线程，执行不同的指令，这叫做线程束的分化。我们已经知道，在一个线程束中所有线程在每个周期中必须执行相同的指令，所以线程束分化似乎会产生一个悖论。

  如果一个线程束中的线程产生分化，线程束将连续执行每一个分支路径，而禁用不执行这一路径的线程。线程束分化会导致性能明显地下降。在前面的例子中可以看到，线程束中并行线程的数量减少了一半：只有16个线程同时活跃地执行，而其他16个被禁用了。条件分支越多，并行性削弱越严重。

  注意，线程束分化只发生在同一个线程束中。在不同的线程束中，不同的条件值不会引起线程束分化。下图显示了线程束分化。

  为了获得最佳的性能，应该避免在同一线程束中有不同的执行路径。请记住，在一个线程块中，线程的线程束分配是确定的。因此，以这样的方式对数据进行分区是可行的（尽管不是微不足道的，但取决于算法），以确保同一个线程束中的所有线程在一个应用程序中使用同一个控制路径。

  现在我们来举个例子，即通过奇偶进入不同的分支执行逻辑，看下如何通过代码来避免线程束分化。先来写一个会造成分化的核函数。我们假设只配置一个x=64的一维线程块，那么只有两个个线程束，线程束内奇数线程（threadIdx.x为奇数）会执行else，偶数线程执行if，分化很严重。

__global__ void mathKernel1(float *c){
    int tid = blockIdx.x * blockDim.x + threadIdx.x;
    float a, b;
    a = b = 0.0f;

    if (tid % 2 == 0){
        a = 100.0f;
    }
    else{
        b = 200.0f;
    }
    c[tid] = a + b;
}

  如果使用线程束方法（而不是线程方法）来交叉存取数据，可以避免线程束分化，并且设备的利用率可达到100%。条件（tid/warpSize）%2==0使分支粒度是线程束大小的倍数；偶数编号的线程执行if子句，奇数编号的线程执行else子句。这个核函数产生相同的输出，但是顺序不同。
第一个线程束内的线程编号tid从0到31，tid/warpSize都等于0，那么就都执行if语句。
第二个线程束内的线程编号tid从32到63，tid/warpSize都等于1，执行else
线程束内没有分支，效率较高。

__global__ void mathKernel2(float *c){
    int tid = blockIdx.x * blockDim.x + threadIdx.x;
    float a, b;
    a = b = 0.0f;

    if ((tid / warpSize) % 2 == 0){
        a = 100.0f;
    }
    else{
        b = 200.0f;
    }
    c[tid] = a + b;
}

完整代码：https://github.com/hujianbin03/dive-into-cuda

#include <cuda_runtime.h>
#include "../include/utils.h"

__global__ void warmingup(float *c){
    int tid = blockIdx.x * blockDim.x + threadIdx.x;
    float a, b;
    a = b = 0.0f;

    if((tid/warpSize) % 2 == 0){
        a = 100.0f;
    }
    else{
        b = 200.0f;
    }
    c[tid] = a + b;
}

__global__ void mathKernel1(float *c){
    int tid = blockIdx.x * blockDim.x + threadIdx.x;
    float a, b;
    a = b = 0.0f;

    if (tid % 2 == 0){
        a = 100.0f;
    }
    else{
        b = 200.0f;
    }
    c[tid] = a + b;
}

__global__ void mathKernel2(float *c){
    int tid = blockIdx.x * blockDim.x + threadIdx.x;
    float a, b;
    a = b = 0.0f;

    if ((tid / warpSize) % 2 == 0){
        a = 100.0f;
    }
    else{
        b = 200.0f;
    }
    c[tid] = a + b;
}

__global__ void mathKernel3(float *c){
    int tid = blockIdx.x * blockDim.x + threadIdx.x;
    float a, b;
    a = b = 0.0f;

    bool ipred = (tid % 2 == 0);
    if (ipred){
        a = 100.0f;
    }
    else{
        b = 200.0f;
    }
    c[tid] = a + b;
}

__global__ void mathKernel4(float *c){
    int tid = blockIdx.x * blockDim.x + threadIdx.x;
    float a, b;
    a = b = 0.0f;

    int itid = tid >> 5;

    if (itid & 0x01 == 0)
    {
        a = 100.0f;
    }
    else
    {
        b = 200.0f;
    }

    c[tid] = a + b;
}

int main(int argc, char **argv){
    // 设置设备
    int dev = 0;
    cudaDeviceProp deviceProp;
    cudaGetDeviceProperties(&deviceProp, dev);
    printf("%s 使用设备 %d: %s\n", argv[0], dev, deviceProp.name);

    // 设置数据
    int size = 64;
    int blocksize = 64;
    if(argc > 1) blocksize = atoi(argv[1]);
    if(argc > 2) size      = atoi(argv[2]);
    printf("数据大小： %d\n", size);

    // 设置线程格、块
    dim3 block (blocksize, 1);
    dim3 grid  ((size+block.x-1)/block.x,1);
    printf("内核配置为：(block %d grid %d)\n", block.x, grid.x);

    // gpu申请内存
    float *d_C;
    size_t nBytes = size * sizeof(float);
    cudaMalloc((float**)&d_C, nBytes);

    // 执行warmup消除开销
    double iStart, iElaps;
    cudaDeviceSynchronize();
    iStart = seconds();
    warmingup<<<grid, block>>>(d_C);
    cudaDeviceSynchronize();
    iElaps = seconds() - iStart;
    printf("warmup          <<< %4d %4d >>> 消耗时间 %lf sec\n", grid.x, block.x, iElaps);

    // 执行mathkernel1
    iStart = seconds();
    mathKernel1<<<grid, block>>>(d_C);
    cudaDeviceSynchronize();
    iElaps = seconds() - iStart;
    printf("mathkernel1     <<< %4d %4d >>> 消耗时间 %lf sec\n", grid.x, block.x, iElaps);

    // 执行mathkernel2
    iStart = seconds();
    mathKernel2<<<grid, block>>>(d_C);
    cudaDeviceSynchronize();
    iElaps = seconds() - iStart;
    printf("mathKernel2     <<< %4d %4d >>> 消耗时间 %lf sec\n", grid.x, block.x, iElaps);

     // 执行mathkernel3
    iStart = seconds();
    mathKernel3<<<grid, block>>>(d_C);
    cudaDeviceSynchronize();
    iElaps = seconds() - iStart;
    printf("mathKernel3     <<< %4d %4d >>> 消耗时间 %lf sec\n", grid.x, block.x, iElaps);

    // 执行mathkernel4
    iStart = seconds();
    mathKernel4<<<grid, block>>>(d_C);
    cudaDeviceSynchronize();
    iElaps = seconds() - iStart;
    printf("mathKernel4     <<< %4d %4d >>> 消耗时间 %lf sec\n", grid.x, block.x, iElaps);

    //释放内存
    cudaFree(d_C);
    cudaDeviceReset();
    return EXIT_SUCCESS;    
}

执行结果如下：

  代码中warmup部分是提前启动一次GPU，因为第一次启动GPU时会比第二次速度慢一些。现在来使用nvprof分析器，可以从GPU中获得指标，从而可以直接观察到线程束分化。

nvprof --metrics branch_efficiency ./simpleDivergence

  分支效率：为未分化的分支与全部分支之比，可以看到这里面所有kernel的分支效率都是100%，可以使用以下公式来计算：：

$$\text{分支效率}=\text{100}\times \frac{\text{分支数 − 分化分支数}}{\text{分支数}}$$

  奇怪的是，没有报告显示出有分支分化（即分支效率是100%）。这个奇怪的现象是CUDA编译器优化导致的结果，它将短的、有条件的代码段的断定指令取代了分支指令（导致分化的实际控制流指令）。但是下面我们用另一种方式，编译器就不会优化了：

__global__ void mathKernel3(float *c){
    int tid = blockIdx.x * blockDim.x + threadIdx.x;
    float a, b;
    a = b = 0.0f;

    bool ipred = (tid % 2 == 0);
    if (ipred){
        ia = 100.0f;
    }
    if (!ipred){
        ib = 200.0f;
    }
    
    c[tid] = a + b;
}

  我们也可以通过编译选项禁用分值预测功能，这样kernel1和kernel3的效率是相近的。如果使用nvprof,会得到下面的结果，没有优化的结果如下：

nvcc -g -G simpleDivergence.cu -o simpleDivergence

  另外，可以用nvprof获得分支和分化分支的事件计数器，如下所示：

nvprof --events branch,divergent_branch ./simpleDivergence

  CUDA的nvcc编译器仍然是在mathKernel1和mathKernel3上执行有限的优化，以保持分支效率在50%以上。注意，mathKernel2不报告分支分化的唯一原因是它的分支粒度是线程束大小的倍数。此外，把mathKernel1中的if…else语句分离为mathKernel3的多个if语句，可以使分化分支的数量翻倍。

重要提示：

• 当一个分化的线程采取不同的代码路径时，会产生线程束分化
• 不同的if-then-else分支会连续执行
• 尝试调整分支粒度以适应线程束大小的倍数，避免线程束分化
• 不同的分化可以执行不同的代码且无须以牺牲性能为代价