任意长度并行前缀和扫描算法《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;
}