关于CUDA向量求和并行规约最后一个warp的展开理解

最新推荐文章于 2024-04-14 13:36:16 发布

粼粼淇

最新推荐文章于 2024-04-14 13:36:16 发布

阅读量745

点赞数 3

分类专栏： CUDA C/C++ 文章标签： CUDA

本文链接：https://blog.csdn.net/lingyunxianhe/article/details/102249270

版权

C/C++ 同时被 2 个专栏收录

19 篇文章 0 订阅

订阅专栏

CUDA

4 篇文章 0 订阅

订阅专栏

CUDA向量求和并行规约有很多内容涉及，在此只说并行规约最后一个warp的展开理解。我参考的是：

Professional CUDA C中的chapter3中的一个例子。Professional CUDA C写的非常好，层次分明，逐渐深入，一个例子从最原始实施到逐步优化完善；

cuda官网上也有并行规约：https://developer.download.nvidia.cn/assets/cuda/files/reduction.pdf，这个很不错，简明扼要。

网上感觉有好多类似的例子，我还是把整个核函数贴出来吧

__global__ void reduceUnrollWarps8 (int *g_idata, int *g_odata, unsigned int n)
{
    // set thread ID
    unsigned int tid = threadIdx.x;
    unsigned int idx = blockIdx.x * blockDim.x * 8 + threadIdx.x;

    // convert global data pointer to the local pointer of this block
    int *idata = g_idata + blockIdx.x * blockDim.x * 8;

    // unrolling 8 
    if (idx + 7 * blockDim.x < n)
    {
        int a1 = g_idata[idx];
        int a2 = g_idata[idx + blockDim.x];
        int a3 = g_idata[idx + 2 * blockDim.x];
        int a4 = g_idata[idx + 3 * blockDim.x];
        int b1 = g_idata[idx + 4 * blockDim.x];
        int b2 = g_idata[idx + 5 * blockDim.x];
        int b3 = g_idata[idx + 6 * blockDim.x];
        int b4 = g_idata[idx + 7 * blockDim.x];
        g_idata[idx] = a1 + a2 + a3 + a4 + b1 + b2 + b3 + b4;
    }

    __syncthreads();

    // in-place reduction in global memory
    for (int stride = blockDim.x / 2; stride > 32; stride >>= 1)
    {
        if (tid < stride)
        {
            idata[tid] += idata[tid + stride];
        }

        // synchronize within threadblock
        __syncthreads();
    }

    // unrolling warp
    //前面的代码不做解释，比较好理解，网上也有一大堆，下面对warp的展开在网上找了好些资料目前没看到一个说的清楚的，好点的解释一下每条语句后为啥不用__syncthreads();了，原因是这些线程属于一个warp,他们隐式同步。要理解这段代码有两个前提：1、cuda的调度单位是warp,不管你当前的线程的数目是否小于32；2、warp中的线程是顺序执行的；;3、warp中的线程是隐式同步的
   
    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];
    }

    // write result for this block to global mem
    if (tid == 0) g_odata[blockIdx.x] = idata[0];
}

对于上述代码，前面的代码不做解释，比较好理解，网上也有一大堆，下面对warp的展开在网上找了好些资料目前没看到一个说的清楚的，好点的解释一下每条语句后为啥不用__syncthreads();了，原因是这些线程属于一个warp,他们隐式同步。

// unrolling warp
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];
    }

要理解这段代码有3个前提：

1、cuda的调度单位是warp,不管你当前的线程的数目是否小于32；

2、warp中的线程是顺序执行的;

3、warp中的线程是隐式同步的；

1 在这段代码中的意义是：当执行如这条语句vmem[tid] += vmem[tid + 32];时是32个线程同时进行的，后面类似

2 在这段代码中的意义是：当执行如这条语句vmem[tid] += vmem[tid + 16];时是32个线程同时进行的，但0-15号线程先进行完毕后（这是理解这段代码的最重要的一点），16-31号线程才进行计算(后面的计算是没有意义的，对求和没有作用，但调度安排是以32为单位的，必须进行呀），因此这条语句执行完毕后，第一次的32个线程的结果累加到前16个线程中。假设不明白2，以为warp中的线程是无序执行的;那么会存在这样的现象，例如执行顺序：vmem[16] += vmem[16 + 16];然后vmem[0] += vmem[0 + 16];此时vmem[0] += vmem[0 + 16]会多加。后面类似

关于warp中的线程是顺序执行的;我进行了一个小测试如下：

__global__ void checkIndex(void)
{
    printf("threadIdx:%d,blockIdx:%d\n", threadIdx.x,blockIdx.x);
}

int main(int argc, char **argv)
{
    // define total data element
    int nElem = 10;
    // define grid and block structure
    dim3 block(4);
    dim3 grid((nElem + block.x - 1) / block.x);
    // check grid and block dimension from host side
    printf("grid.x %d grid.y %d grid.z %d\n", grid.x, grid.y, grid.z);
    printf("block.x %d block.y %d block.z %d\n", block.x, block.y, block.z);
    // check grid and block dimension from device side
    checkIndex<<<grid, block>>>();
    // reset device before you leave
    cudaDeviceReset();

    return(0);
}

第一次运行结果如下：

grid.x 3 grid.y 1 grid.z 1
block.x 4 block.y 1 block.z 1
threadIdx:0,blockIdx:1
threadIdx:1,blockIdx:1
threadIdx:2,blockIdx:1
threadIdx:3,blockIdx:1
threadIdx:0,blockIdx:2
threadIdx:1,blockIdx:2
threadIdx:2,blockIdx:2
threadIdx:3,blockIdx:2
threadIdx:0,blockIdx:0
threadIdx:1,blockIdx:0
threadIdx:2,blockIdx:0
threadIdx:3,blockIdx:0

第二次运行结果如下：

grid.x 3 grid.y 1 grid.z 1
block.x 4 block.y 1 block.z 1
threadIdx:0,blockIdx:1
threadIdx:1,blockIdx:1
threadIdx:2,blockIdx:1
threadIdx:3,blockIdx:1
threadIdx:0,blockIdx:0
threadIdx:1,blockIdx:0
threadIdx:2,blockIdx:0
threadIdx:3,blockIdx:0
threadIdx:0,blockIdx:2
threadIdx:1,blockIdx:2
threadIdx:2,blockIdx:2
threadIdx:3,blockIdx:2

第三次运行结果如下：

grid.x 3 grid.y 1 grid.z 1
block.x 4 block.y 1 block.z 1
threadIdx:0,blockIdx:0
threadIdx:1,blockIdx:0
threadIdx:2,blockIdx:0
threadIdx:3,blockIdx:0
threadIdx:0,blockIdx:2
threadIdx:1,blockIdx:2
threadIdx:2,blockIdx:2
threadIdx:3,blockIdx:2
threadIdx:0,blockIdx:1
threadIdx:1,blockIdx:1
threadIdx:2,blockIdx:1
threadIdx:3,blockIdx:1