cuda - スラストを使用せずにスレッドごとに複数の要素を持つ並列プレフィックス合計

Question

配列の累積合計を見つけるために包括的スキャンを実行しようとしています。ここでハリスムから与えられたアドバイスに従って、私はここで与えられた手順を使用していますが、それらの著者のアドバイスに従って、各スレッドが1つではなく4つの要素を計算してメモリレイテンシをマスクするコードを記述しようとしています.

パフォーマンスが不可欠であり、マルチストリーム機能が必要なため、推力には近づきません。私はCUBを発見したばかりで、それが私の次の努力になるでしょうが、CUDAをよりよく理解するための演習として、マルチブロックソリューションが欲しいです。

以下のコードは、各ブロックに 4 つのデータ要素を割り当てます。各ブロックには 32 の倍数のスレッドが必要です。私のデータには 128 の倍数のスレッドがあるため、この制限は受け入れられます。4*blockDim.x各ブロックには、要素に加えてワープ間で合計する追加の 32 要素用に十分な共有メモリが割り当てられます。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が理解できない新しいエキゾチックなエラーが発生します。

私が間違っている場所を指摘する助けがあれば幸いです