CUDA避免分支分化

线程束中的条件执行可能引起线程束分化，会导致性能变差。通过重新组织数据的获取模式，可以减少或避免线程束分化。

并行归约问题

如果对一个有N个数据的数组求和，串行代码很容易实现

int sum = 0;
for(int i=0; i<N; i++)
    sum += array[i];

如果对大量数据进行并行计算快速求和，可以用以下方法计算

1. 将输入向量划分到更小的数据块中。
2. 用一个线程计算一个数据块的部分和。
3. 对每个数据块的部分和再求和得到最终结果。

并行归约是一种最常见的并行模式，并且是许多并行算法中的一个关键运算。

并行归约中的分化

举个栗子

#include <iostream>
#include <stdio.h>
#include <cuda_runtime.h>
#include <cuda.h>
#include <Windows.h>
#include <random>

__global__ void compute_sum(int* idata, int* odata, const int size) {
    unsigned int tid = threadIdx.x;
    
    int* p = idata + blockIdx.x * blockDim.x;
    if (tid + blockIdx.x * blockDim.x >= size)
        return;

    for (int stride = 1; stride < blockDim.x; stride *= 2) {
        if ((tid % (2 * stride)) == 0) {
            p[tid] += p[tid + stride];
        }
        //栅栏同步所有线程
        __syncthreads();
    }
    if (tid == 0)
        odata[blockIdx.x] = p[0];
}

void print_time(SYSTEMTIME& start, SYSTEMTIME& end) {
    printf("used %d second, %d ms\n", end.wSecond - start.wSecond, end.wMilliseconds - start.wMilliseconds);
}

int main()
{
    using namespace std;
    
    int dev = 0;
    cudaDeviceProp deviceProp;
    cudaGetDeviceProperties(&deviceProp, dev);
    printf("device %d : %s", dev, deviceProp.name);
    cudaSetDevice(dev);

    bool bResult = false;

    int size = 1 << 24;
    printf("    with array size : %d    \n", size);

    int blocksize = 512;

    dim3 block(blocksize, 1);
    dim3 grid((size + block.x - 1) / block.x, 1);
    printf("block size: %d, grid size : %d\n", block.x, grid.x);

    //开辟主机内存
    size_t bytes = size * sizeof(int);
    int* h_idata = (int*)malloc(bytes);
    int* h_odata = (int*)malloc(grid.x * sizeof(int));
    int* tmp = (int*)malloc(bytes);

    //初始化数组
    for (int i = 0; i < size; i++) {
        h_idata[i] = (int)(rand() & 0xFF);
    }
    memcpy(tmp, h_idata, bytes);

    SYSTEMTIME start, end;
    int gpu_sum = 0;

    //开辟GPU内存
    int* d_idata = NULL, * d_odata = NULL;
    cudaMalloc((void**)&d_idata, size * sizeof(int));
    cudaMalloc((void**)&d_odata, grid.x * sizeof(int));

    //CPU时间
    GetSystemTime(&start);
    int cpu_sum = 0;
    for (int i = 0; i < size; i++)
        cpu_sum += h_idata[i];
    GetSystemTime(&end);
    printf("CPU sum:%d ", cpu_sum);
    print_time(start, end);


    //GPU计算
    GetSystemTime(&start);
    cudaMemcpy(d_idata, h_idata, size * sizeof(int), cudaMemcpyHostToDevice);
    
    compute_sum << <grid, block >> > (d_idata, d_odata, size);
    cudaDeviceSynchronize();
    
    cudaMemcpy(h_odata, d_odata, grid.x * sizeof(int), cudaMemcpyDeviceToHost);
    for (int i = 0; i < grid.x; i++) {
        gpu_sum += h_odata[i];
    }
    GetSystemTime(&end);
    printf("GPU sum:%d ", gpu_sum);
    print_time(start, end);

    free(h_odata);
    free(h_idata);
    free(tmp);

    cudaFree(d_idata);
    cudaFree(d_odata);



    return 0;
}

得到的结果为

改善并行归约的分化

在上面的核函数中，有以下表达式

if((tid%(stride*2))==0)

这会令每次迭代中都有新增的一半线程不符合条件，但是这些线程依旧被调度。通过重新组织每个线程的数组索引来强制ID相邻的线程执行求和操作，这样线程束分化就能被归约了。

改进新的核函数

__global__ void new_compute_sum(int* idata, int* odata, const int size) {
    unsigned tid = threadIdx.x;
    unsigned index = threadIdx.x + blockDim.x * blockIdx.x;

    if (index > size) {
        return;
    }

    int* p = idata + blockDim.x * blockIdx.x;

    for (int stride = 1; stride < blockDim.x; stride *= 2) {
        int idx = tid * stride * 2;
        if (idx < blockDim.x) {
            p[idx] += p[idx + stride];
        }

        __syncthreads();
    }

    if (tid == 0)
        odata[blockIdx.x] = p[0];
}

其中，idx = tid*stride*2 为每个线程设置数组访问索引。每轮都有前一半的线程束在计算，而后一半的线程束什么都不做。

在NVIDIA 1660 SUPER上测试，计算会从10ms降低到7ms左右。

由于线程束的大小一般为32，在最后五轮中，在线程数小于32时，分化又会出现。

交错配对的归约

新的核函数

__global__ void new_new_compute_sum(int* idata, int* odata, const int size) {
    unsigned tid = threadIdx.x;
    unsigned index = threadIdx.x + blockDim.x * blockIdx.x;

    if (index > size) {
        return;
    }

    int* p = idata + blockDim.x * blockIdx.x;

    for (int stride = blockDim.x/2; stride>0; stride /= 2) {
        if (tid < stride)
            p[tid] = p[tid] + p[tid + stride];

        __syncthreads();
    }

    if (tid == 0)
        odata[blockIdx.x] = p[0];
}

初始跨度是线程块大小的一半，然后在每次迭代中减少一半。与上一种相比，工作的线程束与线程块没有变化，但是，每个线程在内存中的读写位置不同。性能也差不太多。