0

我正在尝试执行包容性扫描以查找数组的累积和。按照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 )行在它之前的其他逻辑块之前执行。我认为这在经线中是不可能的。

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

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

4

1 回答 1

2

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

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

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

于 2013-07-07T23:46:20.287 回答