【CUDA】Branch Divergence and Unrolling Loop

一、避免分支发散

1.1 并行规约问题

1.2 并行规约中的发散

二、UNrolling Loops

一、避免分支发散

控制流有时依赖于 thread 索引。同一个warp中，一个条件分支可能导致性能很差。通过重新组织数据获取模式可以减少或避免 warp divergence。具体问题查看下面这篇文章：

【CUDA】Warp解析-CSDN博客https://blog.csdn.net/GG_Bruse/article/details/143772619

1.1 并行规约问题

若要计算一个数组N个元素的和，使用CPU编程实现十分容易

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

若数组中的元素非常多，应用并行计算可以大大提高效率。鉴于加法交换律等性质，这个求和过程可以以元素的任意顺序来进行：

将输入数组切割成很多小的块
用 thread 来计算每个块的和
对这些块的结果再求和得最终结果

数组的切割主旨是：用 thread 求数组中按一定规律配对的的两个元素和，然后将所有结果组合成一个新的数组，然后再次求配对两元素和，多次迭代，直到数组中只有一个结果

比较直观的两种实现方式是：

Neighbored pair：每次迭代都是相邻两个元素求和
Interleaved pair：按一定跨度配对两个元素

对于N个元素的数组，该过程需N - 1次求和、 $log_{2}N$ 步。InterleavedPair的跨度是半个数组长度

上述讲的这类问题术语为 reduction problem。Parallel reduction（并行规约）是指迭代减少操作，是并行算法中非常关键的一种操作

1.2 并行规约中的发散

Neighbored pair

在这个kernel中，有两个global memory array，一个用来存放数组所有数据，另一个用来存放部分和。所有 block 独立的执行求和操作

#include <iostream>
#include <cuda_runtime.h>
#include <cuda_runtime_api.h>
using namespace std;// 1 << 28 == 268435456
__global__ 
void reduceNeighbored(int* inputData, int *outputData, unsigned int N) 
{unsigned int threadIndex = threadIdx.x; // 0 - 511unsigned int index = blockIdx.x * blockDim.x + threadIdx.x; // 0 - 524287, 512, 0 - 511 即 0 - 268435455int* iData = inputData + blockIdx.x * blockDim.x; // inputData + (0 - 268434944) 即指向每个block的起始位置if (index >= N) return;// stride = (1 - 511) 1, 2, 4, ..., 500for (int stride = 1; stride < blockDim.x; stride *= 2) {if ((threadIndex % (2 * stride)) == 0) // 每次取threadIndex为偶数的threadiData[threadIndex] += iData[threadIndex + stride];__syncthreads();}if (threadIndex == 0) outputData[blockIdx.x] = iData[0]; // 放入计算出的数据 
}void GPU_ReduceNeighbored(int* dInputData, int* dOutputData, int size, dim3 grid, dim3 block)
{reduceNeighbored<<<grid, block>>>(dInputData, dOutputData, size);
}

因为无法让所有的block同步，所以最后将所有block的结果送回 host 进行串行计算

主函数代码：

void GPU_ReduceNeighbored(int* inputData, int *outputData, int N, dim3 grid, dim3 block);long long Seconds()
{// 获取当前时间点auto now = std::chrono::high_resolution_clock::now();// 将时间点转换为以毫秒为单位的时间间隔auto duration = std::chrono::duration_cast<std::chrono::milliseconds>(now.time_since_epoch());// 获取毫秒数return duration.count();
}int main(int argc, char** argv)
{// set up deviceint device = 0;cudaDeviceProp deviceProp;cudaGetDeviceProperties(&deviceProp, device);printf("%s starting reduction at device %d : %s\n", argv[0], device, deviceProp.name);cudaSetDevice(device);// initialization int size = 1 << 28;printf("Array Size : %d\n", size);int blockSize = 512;dim3 block(blockSize, 1);dim3 grid((size + block.x - 1) / block.x, 1);printf("grid %d , block %d\n", grid.x, block.x);// allocate host memorysize_t bytes = size * sizeof(int);int* hInputData = (int*) malloc(bytes);int* hOutputData = (int*) malloc(grid.x * sizeof(int));int* tmp = (int*) malloc(bytes);for(int i = 0; i < size; ++i) hInputData[i] = (int)(rand() & 0xFF); // 屏蔽最大两个字节memcpy(tmp, hInputData, bytes);// allocate device memoryint* dInputData = nullptr;int* dOutputData = nullptr;cudaMalloc((void **) &dInputData, bytes);cudaMalloc((void **) &dOutputData, grid.x * sizeof(int));int iStart, iElaps;// kernel: reduceNeighboredcudaMemcpy(dInputData, hInputData, bytes, cudaMemcpyHostToDevice);cudaDeviceSynchronize();iStart = Seconds();GPU_ReduceNeighbored(dInputData, dOutputData, size, grid, block);cudaDeviceSynchronize();iElaps = Seconds() - iStart;cudaMemcpy(hOutputData, dOutputData, grid.x * sizeof(int), cudaMemcpyDeviceToHost);long long gpuSum = 0;for (int i = 0; i < grid.x; i++) gpuSum += hOutputData[i];printf("gpu Neighbored elapsed %d ms gpuSum: %lld <<<grid %d block %d>>>\n", iElaps, gpuSum, grid.x, block.x);// free host memoryfree(hInputData);free(hOutputData);// free device memorycudaFree(dInputData);cudaFree(dOutputData);// reset devicecudaDeviceReset();return 0;
}

if ((tid % (2 * stride)) == 0) 该表达式只对偶数ID的线程为true，所以其导致很高的 divergent warps。第一次迭代只有偶数ID的线程执行了指令，但是所有线程都要被调度；第二次迭代，只有四分之一的 thread 是 active 的，但所有 thread 都要被调度。可以重新组织每个线程对应的数组索引来强制ID相邻的thread来处理求和操作

__global__ 
void reduceNeighboredLess(int* inputData, int* outputData, unsigned int N) 
{unsigned int threadIndex = threadIdx.x; // 0 - 511unsigned int index = blockIdx.x * blockDim.x + threadIdx.x; // 0 - 524287, 512, 0 - 511 即 0 - 268435455int* iData = inputData + blockIdx.x * blockDim.x; // inputData + (0 - 268434944) 即指向每个block的起始位置if(index >= N) return;// stride = (1 - 511) 1, 2, 4, ..., 500for (int stride = 1; stride < blockDim.x; stride *= 2) {int idx = 2 * stride * threadIndex; // 2 * (1 - 500) * (0 - 511)if (idx < blockDim.x) iData[idx] += iData[idx + stride];__syncthreads();  }if (threadIndex == 0) outputData[blockIdx.x] = iData[0];
}void GPU_ReduceNeighboredLess(int* dInputData, int* dOutputData, int size, dim3 grid, dim3 block)
{reduceNeighboredLess<<<grid, block>>>(dInputData, dOutputData, size);
}

int index = 2 * stride * tid; 因为步调乘以了2，下面的语句使用block的前半部分thread来执行求和：if (index < blockDim.x)

对于一个有512个 thread 的 block 而言，前八个 warp 执行第一轮 reduction，剩下八个 warp 什么也不干；第二轮，前四个 warp 执行，剩下十二个什么也不干。就不存在divergence了（divergence只发生于同一个warp）。但最后还是会有divergence，因为这个时候需要执行 threads 已经凑不够一个 warp 了

Interleaved pair

Interleaved Pair模式的初始步调是block大小的一半，每个thread处理相隔半个block的两个数据

__global__ 
void reduceInterleaved(int *inputData, int *outputData, unsigned int N) 
{unsigned int threadIndex = threadIdx.x;unsigned int index = blockIdx.x * blockDim.x + threadIdx.x;int *iData = inputData + blockIdx.x * blockDim.x;if(index >= N) return;for (int stride = blockDim.x / 2; stride > 0; stride >>= 1) {if (threadIndex < stride)iData[threadIndex] += iData[threadIndex + stride];__syncthreads();}if (threadIndex == 0) outputData[blockIdx.x] = iData[0];
}void GPU_ReduceInterleaved(int* dInputData, int* dOutputData, int size, dim3 grid, dim3 block)
{reduceInterleaved<<<grid, block>>>(dInputData, dOutputData, size);
}

步调被初始化为block大小的一半：for (int stride = blockDim.x / 2; stride > 0; stride >>= 1)

下面的语句使得第一次迭代时，block的前半部分thread执行相加操作，第二次是前四分之一，以此类推：if (tid < stride)

二、UNrolling Loops

loop unrolling（取消循环）是用来优化循环减少分支的方法，该方法就是把本应在多次 loop 中完成的操作，尽量压缩到一次 loop。循环体展开程度称为loop unrolling factor（循环展开因子），loop unrolling对顺序数组的循环操作性能有很大影响

for (int i = 0; i < 100; ++i) {a[i] = b[i] + c[i];
}

如下重复一次循环体操作，迭代数目将减少一半：

for (int i = 0; i < 100; i += 2) {a[i] = b[i] + c[i];a[i+1] = b[i+1] + c[i+1];
}

从高级语言层面是无法看出性能提升的原因的，需从 low-level instruction 层面去分析，第二段代码循环次数减少了一半，而循环体两句语句的读写操作的执行在CPU上是可以同时执行互相独立的，所以相对第一段，第二段性能要好

Unrolling 在CUDA编程中意义更重。目标是通过减少指令执行消耗，增加更多的独立指令来提高性能。这样就会增加更多的并行操作从而产生更高的指令和内存带宽（bandwidth）。也就提供了更多的 eligible warps 来帮助 hide instruction / memory latency

在前文的 reduceInterleaved 中，每个block处理一部分数据，将这些数据称为 dataBlock

下面的代码是 reduceInterleaved 的修正版本，每个 block 都是以两个 dataBlock 作为源数据进行操作（前文中，每个 block 处理一个 dataBlock）。这是一种循环分区：每个 thread 作用于多个 dataBlock，并且从每个 dataBlock 中取出一个元素处理

__global__ 
void reduceUnrolling(int *inputData, int *outputData, unsigned int N) 
{unsigned int threadIndex = threadIdx.x;unsigned int index = blockIdx.x * blockDim.x * 2 + threadIdx.x;int *iData = inputData + blockIdx.x * blockDim.x * 2;// unrolling 2 data blocksif (index + blockDim.x < N) inputData[index] += inputData[index + blockDim.x];__syncthreads();for (int stride = blockDim.x / 2; stride > 0; stride >>= 1) {if (threadIndex < stride) iData[threadIndex] += iData[threadIndex + stride];__syncthreads();}if (threadIndex == 0) outputData[blockIdx.x] = iData[0];
}  void GPU_ReduceUnrolling(int* dInputData, int* dOutputData, int size, dim3 grid, dim3 block)
{reduceUnrolling<<<grid.x / 2, block>>>(dInputData, dOutputData, size);
}

注意下面的语句，每个 thread 从相邻的 dataBlock中取数据，这一步实际上就是将两个 dataBlock规约成一个：if (idx + blockDim.x < n) g_idata[idx] += g_idata[idx+blockDim.x];

global array index 也要相应的调整，因为，相对之前的版本，同样的数据，只需要原来一半的 thread 就能解决问题。要注意的是，这样做也会降低warp或block的并行性（因为thread减少）：

由于每个block处理两个data block，所以需要调整grid的配置：

reduceUnrolling2<<<grid.x / 2, block>>>(d_idata, d_odata, size);

同一个 thread 中若能有更多的独立的 load/store 操作，会产生更好的性能，因为这样做 memory latency 能够更好的被隐藏

__syncthreads是用来同步block内部thread的。在reduction kernel中，其被用来在每次循环中保证所有 thread 的写 global memory 的操作都已完成，这样才进行下一阶段的计算

当kernel进行到只需要少于或等32个thread（即一个warp）呢？由于使用的SIMT模式，warp内的thread 是有一个隐式的同步过程的。最后几次迭代可以用下面的语句展开：

if (tid < 32) {volatile int *vmem = idata;vmem[tid] += vmem[tid + 32];vmem[tid] += vmem[tid + 16];vmem[tid] += vmem[tid + 8];vmem[tid] += vmem[tid + 4];vmem[tid] += vmem[tid + 2];vmem[tid] += vmem[tid + 1];
}

warp unrolling避免了__syncthreads同步操作，因为这一步本身就没必要

注意volatile修饰符，其告诉编译器每次执行赋值时必须将 vmem[tid] 的值 store 回 global memory。若不这样，编译器或 cache 可能会优化读写global/shared memory。有了这个修饰符，编译器就会认为这个值会被其他 thread 修改，从而使得每次读写都直接去 memory 而不是去 cache 或者 register

__global__ void reduceUnrollWarps8(int *inputData, int *outputData, unsigned int N) 
{unsigned int threadIndex = threadIdx.x;unsigned int index = blockIdx.x * blockDim.x * 8 + threadIdx.x;int *iData = inputData + blockIdx.x * blockDim.x * 8;// unrolling 8if (index + 7*blockDim.x < N) {int a1 = inputData[index];int a2 = inputData[index + blockDim.x];int a3 = inputData[index + 2 * blockDim.x];int a4 = inputData[index + 3 * blockDim.x];int b1 = inputData[index + 4 * blockDim.x];int b2 = inputData[index + 5 * blockDim.x];int b3 = inputData[index + 6 * blockDim.x];int b4 = inputData[index + 7 * blockDim.x];inputData[index] = a1 + a2 + a3 + a4 + b1 + b2 + b3 + b4;}__syncthreads();for (int stride = blockDim.x / 2; stride > 32; stride >>= 1) {if (threadIndex < stride) iData[threadIndex] += iData[threadIndex + stride];__syncthreads();}// unrolling warpif (threadIndex < 32) {volatile int *vmem = iData;vmem[threadIndex] += vmem[threadIndex + 32];vmem[threadIndex] += vmem[threadIndex + 16];vmem[threadIndex] += vmem[threadIndex + 8];vmem[threadIndex] += vmem[threadIndex + 4];vmem[threadIndex] += vmem[threadIndex + 2];vmem[threadIndex] += vmem[threadIndex + 1];}if (threadIndex == 0) outputData[blockIdx.x] = iData[0];
}           void GPU_ReduceUnrolling8(int* dInputData, int* dOutputData, int size, dim3 grid, dim3 block)
{reduceUnrollWarps8<<<grid.x / 8, block>>>(dInputData, dOutputData, size);
}

若在编译时已知了迭代次数，就可以完全把循环展开。Fermi 和 Kepler 每个 block 的最大 thread 数目都是1024，文中的 kernel 的迭代次数都是基于 blockDim 的，所以完全展开循环是可行的

__global__ 
void reduceCompleteUnrollWarps8(int *inputData, int *outputData, unsigned int N) 
{unsigned int threadIndex = threadIdx.x;unsigned int index = blockIdx.x * blockDim.x * 8 + threadIdx.x;int *iData = inputData + blockIdx.x * blockDim.x * 8;// unrolling 8if (index + 7 * blockDim.x < N) {int a1 = inputData[index];int a2 = inputData[index + blockDim.x];int a3 = inputData[index + 2 * blockDim.x];int a4 = inputData[index + 3 * blockDim.x];int b1 = inputData[index + 4 * blockDim.x];int b2 = inputData[index + 5 * blockDim.x];int b3 = inputData[index + 6 * blockDim.x];int b4 = inputData[index + 7 * blockDim.x];inputData[index] = a1 + a2 + a3 + a4 + b1 + b2 + b3 + b4;}__syncthreads();if (blockDim.x>=1024 && threadIndex < 512) iData[threadIndex] += iData[threadIndex + 512];__syncthreads();        if (blockDim.x>=512 && threadIndex < 256) iData[threadIndex] += iData[threadIndex + 256];__syncthreads();if (blockDim.x>=256 && threadIndex < 128) iData[threadIndex] += iData[threadIndex + 128];__syncthreads();if (blockDim.x>=128 && threadIndex < 64) iData[threadIndex] += iData[threadIndex + 64];__syncthreads();// unrolling warpif (threadIndex < 32) {volatile int *vsmem = iData;vsmem[threadIndex] += vsmem[threadIndex + 32];vsmem[threadIndex] += vsmem[threadIndex + 16];vsmem[threadIndex] += vsmem[threadIndex + 8];vsmem[threadIndex] += vsmem[threadIndex + 4];vsmem[threadIndex] += vsmem[threadIndex + 2];vsmem[threadIndex] += vsmem[threadIndex + 1];}if (threadIndex == 0) outputData[blockIdx.x] = iData[0];
}                void GPU_ReduceCompleteUnrollWarps8(int* dInputData, int* dOutputData, int size, dim3 grid, dim3 block)
{reduceCompleteUnrollWarps8<<<grid.x / 8, block>>>(dInputData, dOutputData, size);
}