在不使用推力的情况下，每个线程具有多个元素的并行前缀总和

Posted 2023-02-16

【中文标题】在不使用推力的情况下，每个线程具有多个元素的并行前缀总和

【英文标题】：Parallel prefix sum with multiple elements per thread without using thrust

【发布时间】：2013-07-07 16:42:30

【问题描述】：

我正在尝试执行包容性扫描以查找数组的累积和。遵循harrism here 给出的建议，我正在使用here 给出的过程，但是按照这些作者的建议，我正在尝试编写代码，让每个线程计算4 个元素而不是一个元素来掩盖内存延迟.

我远离推力，因为性能至关重要，而且我需要多流能力。我刚刚发现了 CUB，这将是我的下一个努力，但我想要一个多块解决方案，并且还想知道我在现有代码上哪里出了问题，作为更好地理解 CUDA 的练习。

下面的代码为每个块分配 4 个数据元素，其中每个块必须有 32 个线程的倍数。我的数据将有 128 个线程的倍数，所以这个限制对我来说是可以接受的。为4*blockDim.x 元素的每个块分配了足够的共享内存，再加上额外的 32 个元素以在 warp 之间求和。 scanBlockAnyLength 然后添加必要的偏移量以纠正扭曲之间的不匹配，将每个扭曲的最终值保存到设备全局内存中的dev_blockSum。 sumWarp4_32 然后扫描这个数组找到最终纠正块之间的不匹配，然后添加到kernel_sumBlock

#include<cuda.h>
#include<iostream>
using std::cout;
using std::endl;

#define MAX_THREADS 1024
#define MAX_BLOCKS 65536
#define N 512

__device__ float sumWarp4_128(float* ptr, const int tidx = threadIdx.x) 
    const unsigned int lane = tidx & 31;
    const unsigned int warpid = tidx >> 5; //32 threads per warp

    unsigned int i = warpid*128+lane; //first element of block data set this thread looks at

    if( lane >= 1 ) ptr[i] += ptr[i-1];
    if( lane >= 2 ) ptr[i] += ptr[i-2];
    if( lane >= 4 ) ptr[i] += ptr[i-4];
    if( lane >= 8 ) ptr[i] += ptr[i-8];
    if( lane >= 16 ) ptr[i] += ptr[i-16];

    if( lane==0 ) ptr[i+32] += ptr[i+31];

    if( lane >= 1 ) ptr[i+32] += ptr[i+32-1];
    if( lane >= 2 ) ptr[i+32] += ptr[i+32-2];
    if( lane >= 4 ) ptr[i+32] += ptr[i+32-4];
    if( lane >= 8 ) ptr[i+32] += ptr[i+32-8];
    if( lane >= 16 ) ptr[i+32] += ptr[i+32-16];

    if( lane==0 ) ptr[i+64] += ptr[i+63];

    if( lane >= 1 ) ptr[i+64] += ptr[i+64-1];
    if( lane >= 2 ) ptr[i+64] += ptr[i+64-2];
    if( lane >= 4 ) ptr[i+64] += ptr[i+64-4];
    if( lane >= 8 ) ptr[i+64] += ptr[i+64-8];
    if( lane >= 16 ) ptr[i+64] += ptr[i+64-16];

    if( lane==0 ) ptr[i+96] += ptr[i+95];

    if( lane >= 1 ) ptr[i+96] += ptr[i+96-1];
    if( lane >= 2 ) ptr[i+96] += ptr[i+96-2];
    if( lane >= 4 ) ptr[i+96] += ptr[i+96-4];
    if( lane >= 8 ) ptr[i+96] += ptr[i+96-8];
    if( lane >= 16 ) ptr[i+96] += ptr[i+96-16];

    return ptr[i+96];

__host__ __device__ float sumWarp4_32(float* ptr, const int tidx = threadIdx.x) 
    const unsigned int lane = tidx & 31;
    const unsigned int warpid = tidx >> 5; //32 elements per warp

    unsigned int i = warpid*32+lane; //first element of block data set this thread looks at

    if( lane >= 1 ) ptr[i] += ptr[i-1];
    if( lane >= 2 ) ptr[i] += ptr[i-2];
    if( lane >= 4 ) ptr[i] += ptr[i-4];
    if( lane >= 8 ) ptr[i] += ptr[i-8];
    if( lane >= 16 ) ptr[i] += ptr[i-16];

    return ptr[i];

__device__ float sumBlock4(float* ptr, const int tidx = threadIdx.x, const int bdimx = blockDim.x ) 
    const unsigned int lane = tidx & 31;
    const unsigned int warpid = tidx >> 5; //32 threads per warp

    float val = sumWarp4_128(ptr);
    __syncthreads();//should be included

    if( tidx==bdimx-1 ) ptr[4*bdimx+warpid] = val;
    __syncthreads();

    if( warpid==0 ) sumWarp4_32((float*)&ptr[4*bdimx]);
    __syncthreads();

    if( warpid>0 ) 
        ptr[warpid*128+lane] += ptr[4*bdimx+warpid-1];
        ptr[warpid*128+lane+32] += ptr[4*bdimx+warpid-1];
        ptr[warpid*128+lane+64] += ptr[4*bdimx+warpid-1];
        ptr[warpid*128+lane+96] += ptr[4*bdimx+warpid-1];
    
    __syncthreads();
    return ptr[warpid*128+lane+96];

__device__ void scanBlockAnyLength4(float *ptr, float* dev_blockSum, const float* dev_input, float* dev_output, const int idx = threadIdx.x, const int bdimx = blockDim.x, const int bidx = blockIdx.x) 

    const unsigned int lane = idx & 31;
    const unsigned int warpid = idx >> 5;

    ptr[lane+warpid*128] = dev_input[lane+warpid*128+bdimx*bidx*4];
    ptr[lane+warpid*128+32] = dev_input[lane+warpid*128+bdimx*bidx*4+32];
    ptr[lane+warpid*128+64] = dev_input[lane+warpid*128+bdimx*bidx*4+64];
    ptr[lane+warpid*128+96] = dev_input[lane+warpid*128+bdimx*bidx*4+96];
    __syncthreads();

    float val = sumBlock4(ptr);
    __syncthreads();
    dev_blockSum[0] = 0.0f;
    if( idx==0 ) dev_blockSum[bidx+1] = ptr[bdimx*4-1];

    dev_output[lane+warpid*128+bdimx*bidx*4] = ptr[lane+warpid*128];
    dev_output[lane+warpid*128+bdimx*bidx*4+32] = ptr[lane+warpid*128+32];
    dev_output[lane+warpid*128+bdimx*bidx*4+64] = ptr[lane+warpid*128+64];
    dev_output[lane+warpid*128+bdimx*bidx*4+96] = ptr[lane+warpid*128+96];
    __syncthreads();

__global__ void kernel_sumBlock(float* dev_blockSum, const float* dev_input, float*   dev_output ) 
    extern __shared__ float ptr[];
    scanBlockAnyLength4(ptr,dev_blockSum,dev_input,dev_output);

__global__ void kernel_offsetBlocks(float* dev_blockSum, float* dev_arr) 
    const int tidx = threadIdx.x;
    const int bidx = blockIdx.x;
    const int bdimx = blockDim.x;

    const int lane = tidx & 31;
    const int warpid = tidx >> 5;
    if( warpid==0 ) sumWarp4_32(dev_blockSum);
    float val = dev_blockSum[warpid];
    dev_arr[warpid*128+lane] += val;
    dev_arr[warpid*128+lane+32] += val;
    dev_arr[warpid*128+lane+64] += val;
    dev_arr[warpid*128+lane+96] += val;

void scan4( const float input[], float output[]) 
    int blocks = 2;
    int threadsPerBlock = 64; //multiple of 32
    int smemsize = (threadsPerBlock*4+32)*sizeof(float);

    float* dev_input, *dev_output;
    cudaMalloc((void**)&dev_input,blocks*threadsPerBlock*4*sizeof(float));
    cudaMalloc((void**)&dev_output,blocks*threadsPerBlock*4*sizeof(float));

    float *dev_blockSum;
    cudaMalloc((void**)&dev_blockSum,blocks*sizeof(float));

    int offset = 0;
    int Nrem = N;
    int chunksize;
    while( Nrem ) 
        chunksize = max(Nrem,blocks*threadsPerBlock*4);
        cudaMemcpy(dev_input,(void**)&input[offset],chunksize*sizeof(float),cudaMemcpyHostToDevice);
        kernel_sumBlock<<<blocks,threadsPerBlock,smemsize>>>(dev_blockSum,dev_input,dev_output);
        kernel_offsetBlocks<<<blocks,threadsPerBlock>>>(dev_blockSum,dev_output);
        cudaMemcpy((void**)&output[offset],dev_output,chunksize*sizeof(float),cudaMemcpyDeviceToHost);
        offset += chunksize;
        Nrem -= chunksize;
    
    cudaFree(dev_input);
    cudaFree(dev_output);


int main() 
    float h_vec[N], sol[N];
    for( int i = 0; i < N; i++ ) h_vec[i] = (float)i+1.0f;

    scan4(h_vec,sol);

    cout << "solution:" << endl;
    for( int i = 0; i < N; i++ ) cout << i << " " << (i+2)*(i+1)/2 << " " << sol[i] << endl;
    return 0;

在我看来，代码会抛出错误，因为 sumWarp4_128 中的行在扭曲中没有按顺序执行。即，if( lane==0 ) 行在它之前的其他逻辑块之前执行。我认为这在经线中是不可能的。

如果我在lane==0 调用之前和之后__syncthreads()，我会收到一些我无法弄清楚的新奇错误。

任何帮助指出我出错的地方将不胜感激

【问题讨论】：

如果您使用cuda-memcheck 运行代码，您会发现报告了几种不同类型的错误。特别是，您在内核sumBlock 中有越界访问。我将首先在所有 cuda API 调用和内核调用上正确实现 cuda error checking，然后整理出正在生成的错误。如果您需要有关细节方面的帮助，请发布您收到的错误输出以及您到目前为止所做的工作以解决问题。

【参考方案1】：

由于共享数据的线程之间没有同步，您正在编写的代码存在竞争条件。虽然这确实可以在当前硬件上完成，以便在 warp 内进行通信（所谓的 warp 同步编程），但非常不鼓励这样做，因为代码中的竞争条件可能会导致它在未来可能的硬件上失败。

虽然通过每个线程处理多个项目确实可以获得更高的性能，但 4 并不是一个神奇的数字——如果可能，您应该将其设为可调参数。例如，CUDPP 每个线程使用 8 个。

我强烈建议您为此使用 CUB。您应该使用cub::BlockLoad() 来为每个线程加载多个项目，并使用cub::BlockScan() 来扫描它们。然后你只需要一些代码来组合多个块。最节省带宽的方法是使用 Thrust 使用的“Reduce-Scan-Scan”方法。首先减少每个块（cub::BlockReduce）并将每个块的总和存储到blockSums 数组中。然后扫描该数组以获取每个块的偏移量。然后对块执行 cub::BlockScan 并将先前计算的每块偏移添加到每个元素。

【讨论】：

FWIW，here's a reference implementation of scan。它采用了 Mark 提到的 reduce-scan-scan 方法。它不使用 CUB，但速度很快。

它的模板非常深入，很难找到核心功能的实现位置。最低级别的扫描功能在哪里定义？您指的是哪个版本“快”——sean_scan 或 my_scan？

我对 CUB 的存在一无所知。它是否使用多个流进行数据传输和计算？我正在处理的实际问题使用的是自引用的多 GB 数据集，这就是为什么我认为 Thrust 不适合我的工具。快速检查 CUB 主页并不清楚编辑：每次扫描都在一个易于放入设备内存的数组上，但计算要扫描的数组需要从主机获取其他数组

CUB 是一种设备端 API，它与流和主机设备数据传输正交且不关心。