任意长度并行前缀和 扫描算法 《PMPP》笔记

下面的算法针对于任意长度输入

对于大数据集，首先将输入分为几段，每一段放进共享内存并用一个线程块处理，比如一个线程块使用1024个线程的话，每个块最多能处理2048个元素。

在前面代码中，一个块最后的执行结果保存到了Y数组中，Y 数组保存了每个段扫描的结果，可以称之为扫描块， 一个扫描块只保存了当前块中前面所有元素的累加值，需要把这些扫描块合并到一个最终的结果中。

上述栗子， 在16个输入的数组中，分为4个扫描块，kernel将4个扫描块看做独立的输入数据集处理，扫描kernel结束之后，每个Y元素保存了这个扫描块中扫描的结果。

每个扫描块最后一个元素时当前扫描块中输入元素的总和。

在第二步中，从每个扫描块中收集最后一个元素，放进一个数组S中，然后对此数组进行扫描，然后将扫描S数组后的值累加到对应的扫描块上。

可以使用3个kernel实现层级扫描，第一个kernel和之前的kernel没有太大差别（都是针对块内进行扫描）， 需要添加一个中间变量S，其维度为 inputSize/SECTION_SIZE, 在kernel的最后，需要块的最后一个线程把当前扫描块中最后值写到S中blockIdx.x 位置上。

第二个kernel和之前的kernel也一样，只是使用S作为输入，修改S的内容并将之作为输出。

第三个kernel接受S和Y数组作为输入，然后将输出写回到Y， 将一个S的元素加到对应扫描块的Y元素上。

/*
处理任意长度输入的并行归约， 包括3个层级kernel
*/
__global__ void tier1_scan_kernel(
    float* dev_x, float *dev_y, float *dev_s, unsigned int inputSize){
    // 第一层级，实现每个块内的归约，并将归约后的最后一个元素写到S中
    __shared__ float XY[SECTION_SIZE];
    int idx = blockIdx.x * blockDim.x +threadIdx.x;
    if(idx < inputSize){
        XY[threadIdx.x] = dev_x[idx];
    }

    // 归约阶段
    for(unsigned int stride=1;stride<blockDim.x; stride*=2){
        __syncthreads();
        int index = (threadIdx.x+1)*2*stride - 1;
        if(index<blockDim.x){
            XY[index] += XY[index-stride];
        }
    }

    // 分发阶段
    for(int stride=SECTION_SIZE/4; stride>0; stride/=2){
        __syncthreads();
        int index = (threadIdx.x+1)*stride*2 - 1;
        if(index+stride< SECTION_SIZE){
            XY[index+stride] += XY[index];
        }
    }

    __syncthreads();
    dev_y[idx] = XY[threadIdx.x];
    if (threadIdx.x == 0){
        dev_s[blockIdx.x] = XY[SECTION_SIZE-1];
    }
}

__global__ void tier2_scan_kernel(float * dev_s, unsigned int inputSize){
    __shared__ float XY[SECTION_SIZE];
    int idx = blockIdx.x * blockDim.x +threadIdx.x;
    if(idx < inputSize){
        XY[threadIdx.x] = dev_s[idx];
    }

    // 归约阶段
    for(unsigned int stride=1;stride<blockDim.x; stride*=2){
        __syncthreads();
        int index = (threadIdx.x+1)*2*stride - 1;
        if(index<blockDim.x){
            XY[index] += XY[index-stride];
        }
    }

    // 分发阶段
    for(int stride=SECTION_SIZE/4; stride>0; stride/=2){
        __syncthreads();
        int index = (threadIdx.x+1)*stride*2 - 1;
        if(index+stride< SECTION_SIZE){
            XY[index+stride] += XY[index];
        }
    }

    __syncthreads();
    dev_s[idx] = XY[threadIdx.x];
}

__global__ void tier3_scan_kernel(float *dev_y, float *dev_s, unsigned int inputSize){
    int idx = blockIdx.x * blockDim.x + threadIdx.x;
    if (idx< inputSize){
        dev_y[idx] += dev_s[blockIdx.x];
    }
}

void func_scan_gpu3(float* x, unsigned int length){
    float *y = new float[length];

    float *dev_x, *dev_y, *dev_s;
    cudaMalloc((void**)&dev_x, length*sizeof(float));
    cudaMalloc((void**)&dev_y, length*sizeof(float));
    unsigned int blocks = (length + SECTION_SIZE -1)/ SECTION_SIZE;

    cudaMemcpy(dev_x, x, length*sizeof(float), cudaMemcpyHostToDevice);
    cudaMalloc((void**)&dev_s, blocks*sizeof(float));
    tier1_scan_kernel<<<blocks, SECTION_SIZE>>>(dev_x, dev_y, dev_s, length);
    tier2_scan_kernel<<<1, blocks>>>(dev_s, blocks);
    tier3_scan_kernel<<<blocks, SECTION_SIZE>>>(dev_y,dev_s, length);

    cudaMemcpy(y, dev_y,length*sizeof(float), cudaMemcpyDeviceToHost);
    print1DArr(y, SECTION_SIZE);

    cudaFree(dev_x);
    cudaFree(dev_y);
    cudaFree(dev_s);
    delete[] y;
}