GPU实现Prefix Sum Scan (cumulative sum )算法

weixin_42849849

已于 2022-01-28 20:12:44 修改

阅读量787

点赞数

文章标签：算法 CUDA GPU

于 2021-04-12 19:07:46 首次发布

本文链接：https://blog.csdn.net/weixin_42849849/article/details/115634494

版权

背景

对给定NBATCH组数据

float Array[NBATCH][LENGTH];

对每组数据需要进行一次forward正向scan sum运算，然后进行backward反向一次scan sum。这种操作来源于地震道数据在偏移前的反假频滤波需要。其中NBATCH数量级可在10k左右，LENGTH数量级在10k左右。

现有Thrust或CUB库没有提供适合多个数组的批量操作，对正反向结合也需要较复杂的逻辑表达。此外对于LENGTH较大的数组，需要考虑单精度数累加过程中精度丢失问题(rounding error)，在内部运算中采用两种方式：一种是用double作为临时累加求和计算；另一种用float2模拟double的求和运算，以期避免精度丢失问题。

技术路线

Thrust或CUB。针对本问题，构建有难度
基于shared memory优化版: gpu_scansum_smem_v0。聚合访问global memory, 设置PAD避免bank conflict
基于shared memory优化版: gpu_scansum_smem_v1。移除或减少分支，避免Warp内同步
基于Warp Shuffles(__shfl_up_sync)原语: gpu_scansum_warp_v0。没有使用shared memory,直接在warp内线程的寄存器进行数据交换和计算
针对float,double,float2三个累加求和类型，对kernel进行模板化，以便于性能测试。 float2模拟double计算基于Lewis的书内源码，进行少量改装以适合C++模板!
编译时不能打开浮点计算优化选项，保证计算满足IEEE754标准！

测试环境

环境#1

系统: Ubuntu18.04
GCC: 7.5
CUDA Toolkit: 11.2
GPU: GeForce GTX 1070 Ti (计算能力6.1, Pascal架构)
优化选项: -O3

环境#2

系统: win10
vs2019
CUDA Toolkit: 11.1
GPU: GeForce GTX TITAN X(计算能力5.2，Maxwell架构)

实现

cpu上

    template<class T>
    void  cpu_scansum_kernel(float *h_data, int  len, int nbatch)
    {
        for(int i=0; i<nbatch; i++)
        {
            T sum=0.0f;
            //forward
            for(int j=0; j<len; j++)
            {
                sum=sum+(h_data+i*len)[j];
                (h_data+i*len)[j]=(float)sum;
            }

            sum=0.0f;
            //backward
            for(int j=len-1;j>=0; j--)
            {
                sum=sum+(h_data+i*len)[j];
                (h_data+i*len)[j]=(float)sum;
            }
        }
    }

GPU: gpu_scansum_smem_v0。此实现一个warp处理32个数组，每次将每个数组的32个数加载到shared memory,然后在shared memory内每个thread进行累加求和。这里要考虑数组数目和长度不是32整数倍情况，所以内部会有边界判断。shared memory考虑到bank conflict情况，在第二维度上设为32+PAD, PAD=1时候可以避免bank conflict, PAD=0时候bank conflict最严重!


    template<typename T>
    __global__ void  gpu_cumsum_smem_v0(float *d_data, int len, int  nbatch)
    {
        //threads(32,1,1):
        //assert(blockDim.x==32&&blockDim.y==1 && blockDim.z==1)
	    
	    static const int PAD=1;
        __shared__ float my_floats[32][32+PAD];

        const int tid=threadIdx.x;
        const int bid=blockIdx.x;
        const int  blocks=gridDim.x;
        for(int i=32*bid; i<nbatch; i+=32*blocks)
        {   
            //forward scan sum (cumsum)
            T sum=0.0f;
            for(int j=0; j<len; j+=32)
            {
                //load from global mem to shared mem
                if(tid+j<len)
                {
                    for(int k=0; k<32; k++)
                    {
                        if(i+k<nbatch)
                            my_floats[k][tid]=(d_data+(i+k)*len)[j+tid];
                    }
                }
                //__syncthreads();
               for(int k=0; k<32; k++)
               {
                   sum=sum+my_floats[tid][k];
                   my_floats[tid][k]=(float)sum;
               }
                
                //__syncthreads();
                
                //store shared mem to global mem
                if(tid+j<len)
                {
                    //for(uint32_t k=0; k<32; k++)
                    for(int k=0; k<32; k++)
                    {
                        if(i+k<nbatch)
                            (d_data+(i+k)*len)[j+tid]=my_floats[k][tid];
                    }
                }
            }

            //backward scan sum (cumsum)
            sum=0.0f;
            for(int j=len-1; j>=0; j-=32)
            {
                //load from global mem to shared mem
                if(j>=tid)
                {
                    for(int k=0; k<32; k++)
                    {
                        if(i+k<nbatch)
                            my_floats[k][tid]=(d_data+(i+k)*len)[j-tid];
                    }
                }
                //__syncthreads();
                //scan sum op
                for(int k=0; k<32; k++)
                {
                    sum=sum+my_floats[tid][k];
                    my_floats[tid][k]=(float)sum;
                }

                //__syncthreads();
                
                //store shared mem to global mem
                if(j>=tid)
                {
                    for(int k=0; k<32;k++)
                    {
                        if(i+k<nbatch)
                            (d_data+(i+k)*len)[j-tid]=my_floats[k][tid];
                    }
                }
            }
        }
    }

GPU: gpu_scansum_smem_v1。相比较gpu_scansum_smem_v0版本，考虑到nbatch和len不是32整数倍情况，将循环内部边界判断if分支剥离出来，单独处理，以期望提高性能。实际测量对性能提高很小，可能和此问题是非计算密集型运算有关！对循环边界处理使代码明显变长

    template<typename T>
    __global__ void  gpu_cumsum_smem_v1(float *d_data, int len, int  nbatch)
    {
        __shared__ float my_floats[32][32+PAD];

        const int tid=threadIdx.x;
        const int bid=blockIdx.x;
        const int blocks=gridDim.x;

        int nbatch32=nbatch>>5;
        int nbatch_left=nbatch&31;
        int len32=len&~31;
        int len_left=len&31;
        int tid_lt_len_left=tid<len_left;
        int tid_lt_nbatch_left=tid<nbatch_left;

#define FORWARD(j,tid) (j+tid)
#define BACKWARD(j,tid) (j-tid)
       
        int i;
        for(i=bid; i<nbatch32; i+=blocks)
        {
            float *array=i*32*len+d_data;
            /********************************************
            *  forward scan sum
            ***********************************************/
            T sum=0.0f;
            int j;
            for(j=0; j<len32; j+=32)
            {
                #pragma unroll 
                for(int k=0; k<32; k++)
                {
                    my_floats[k][tid]=(array+k*len)[FORWARD(j,tid)];
                }          
                #pragma unroll 
                for(int k=0; k<32; k++)
                {
                    sum=sum+my_floats[tid][k];
                    my_floats[tid][k]=(float)sum;
                }
                #pragma unroll 
                for(int k=0; k<32; k++)
                {
                    (array+k*len)[FORWARD(j,tid)]=my_floats[k][tid];
                }
            }
            
            if(len_left==0) goto STEP_BACKWARD;

            // handle len_left>0
            //if(tid<len_left)
            if(tid_lt_len_left)
            {
                #pragma unroll 
                for(int k=0; k<32; k++)
                {
                    my_floats[k][tid]=(array+k*len)[FORWARD(j,tid)];
                }
            }

            for(int k=0; k<len_left; k++)
            {
                sum=sum+my_floats[tid][k];
                my_floats[tid][k]=(float)sum;
            }

            //if(tid<len_left)
            if(tid_lt_len_left)
            {
                #pragma unroll 
                for(int k=0; k<32; k++)
                    (array+k*len)[FORWARD(j,tid)]=my_floats[k][tid];

            }

            STEP_BACKWARD:

            /***********************************************************
            *  backward scan sum
            ************************************************************/

            sum=0.0f;
            for(j=len-1; j>=len_left; j-=32)
            {
                #pragma unroll 
                for(int k=0; k<32; k++)
                {
                    my_floats[k][tid]=(array+k*len)[BACKWARD(j,tid)];
                }
                #pragma unroll 
                for(int k=0; k<32; k++)
                {
                    sum=sum+my_floats[tid][k];
                    my_floats[tid][k]=(float)sum;
                }
                #pragma unroll 
                for(int k=0; k<32; k++)
                {
                    (array+k*len)[BACKWARD(j,tid)]=my_floats[k][tid];
                }
            }

            if(len_left==0) continue;

            // handle len_left>0
            //if(tid<len_left)
            if(tid_lt_len_left)
            {
                #pragma unroll 
                for(int k=0; k<32; k++)
                {
                    my_floats[k][tid]=(array+k*len)[BACKWARD(j,tid)];
                }
            }

            for(int k=0; k<len_left; k++)
            {
                sum=sum+my_floats[tid][k];
                my_floats[tid][k]=(float)sum;
            }

            //if(tid<len_left)
            if(tid_lt_len_left)
            {
                #pragma unroll 
                for(int k=0; k<32; k++)
                    (array+k*len)[BACKWARD(j,tid)]=my_floats[k][tid];
            }

        }

        //handle nbatch_left>0
        //if(bid==0 && nbatch_left)
        if(bid==blocks-1 && nbatch_left>0)
        {
            float *array=d_data+32*nbatch32*len;
            T sum=0.0f; 
            int j;
            for(j=0; j<len32; j+=32)
            {
                for(int k=0; k<nbatch_left; k++)
                {
                    my_floats[k][tid]=(array+k*len)[FORWARD(j,tid)];
                }

                //if(tid<nbatch_left)
                if(tid_lt_nbatch_left)
                {
                    #pragma unroll 
                    for(int k=0; k<32; k++)
                    {
                        sum=sum+my_floats[tid][k];
                        my_floats[tid][k]=(float)sum;
                    }
                }

                for(int k=0; k<nbatch_left; k++)
                {
                    (array+k*len)[FORWARD(j,tid)]=my_floats[k][tid];
                }
            }

            //handle len_left>0

            //if(tid<len_left)
            if(tid_lt_len_left)
            {
                for(int k=0; k<nbatch_left; k++)
                {
                    my_floats[k][tid]=(array+k*len)[FORWARD(j,tid)];
                }
            }

            //if(tid<nbatch_left)
            if(tid_lt_nbatch_left)
            {
                for(int k=0; k<len_left; k++)
                {
                    sum=sum+my_floats[tid][k];
                    my_floats[tid][k]=(float)sum;
                }
            }
            //if(tid<len_left)
            if(tid_lt_len_left)
            {
                for(int k=0; k<nbatch_left; k++)
                {
                    (array+k*len)[FORWARD(j,tid)]=my_floats[k][tid];
                }
            }

            //backward
            sum=0.0f;
            for(j=len-1; j>=len_left; j-=32)
            {
                for(int k=0; k<nbatch_left; k++)
                {
                    my_floats[k][tid]=(array+k*len)[BACKWARD(j,tid)];
                }

                //if(tid<nbatch_left)
                if(tid_lt_nbatch_left)
                {
                    #pragma unroll 
                    for(int k=0; k<32; k++)
                    {
                        sum=sum+my_floats[tid][k];
                        my_floats[tid][k]=(float)sum;
                    }
                }

                for(int k=0; k<nbatch_left; k++)
                {
                    (array+k*len)[BACKWARD(j,tid)]=my_floats[k][tid];
                }
            }

            //handle len_left>0
            //if(tid<len_left)
            if(tid_lt_len_left)
            {
                for(int k=0; k<nbatch_left; k++)
                {
                    my_floats[k][tid]=(array+k*len)[BACKWARD(j,tid)];
                }
            }

            //if(tid<nbatch_left)
            if(tid_lt_nbatch_left)
            {
                for(int k=0; k<len_left; k++)
                {
                    sum=sum+my_floats[tid][k];
                    my_floats[tid][k]=(float)sum;
                }
            }
            //if(tid<len_left)
            if(tid_lt_len_left)
            {
                for(int k=0; k<nbatch_left; k++)
                {
                    (array+k*len)[BACKWARD(j,tid)]=my_floats[k][tid];
                }
            }            
        }
#undef FORWARD
#undef BACKWARD        
    }

GPU: gpu_scansum_warp_v0。用__shfl_up_sync() 指令原语进行warp内进程寄存器数据交换，实现scan sum计算! __shfl_up_sync对4Byte数据性能可以，对8Byte宽类型延迟明显增加。

    template<typename T>
    __device__  void warpscan_forward(int tid, float *d_data, int n)
    {
        static_assert(sizeof(T)==4 || sizeof(T)==8);
        static_assert(std::is_same<T,float>::value || std::is_same<T,double>::value || std::is_same<T,FLOAT2>::value);

        T pre_sum=0.0f;
        T sum;
        int j=0; 
        for(j=0; j<n-31; j+=32)
        {
            sum=d_data[j+tid];
            for(int d=1; d<32; d*=2)
            {   
                T temp;
   
                if(std::is_same<T,FLOAT2>::value)
                    *(uint64_t*)&temp=__shfl_up_sync(-1,*(uint64_t*)&sum,d);
                else 
                    temp=__shfl_up_sync(-1,sum,d);
   
                if(tid>=d)
                    sum=sum+temp;
            }
            d_data[j+tid]=(float)(sum+pre_sum);

            if(std::is_same<T,FLOAT2>::value)
                *(uint64_t*)&sum=__shfl_sync(-1,*(uint64_t*)&sum,31);
            else 
                sum=__shfl_sync(-1,sum,31);
   
            pre_sum=pre_sum+sum;
        }

        if(j+tid<n)
            sum=d_data[j+tid];

        for(int d=1; d<n-j; d*=2)
        {   
            T temp;

            if(std::is_same<T,FLOAT2>::value)
                *(uint64_t*)&temp=__shfl_up_sync(-1,*(uint64_t*)&sum,d);
            else 
                temp=__shfl_up_sync(-1,sum,d);

            if(tid>=d)
                sum=sum+temp;
        }
        if(j+tid<n)
            d_data[j+tid]=(float)(sum+pre_sum);
   }//warpscan_forward

    template<typename T>
    __device__  void warpscan_backward(int tid, float *d_data, int n)
    {

        static_assert(sizeof(T)==4 || sizeof(T)==8);
        static_assert(std::is_same<T,float>::value || std::is_same<T,double>::value || std::is_same<T,FLOAT2>::value);

        T pre_sum=0.0f;
        T sum;
        int j; 
        for(j=n-1; j>=31; j-=32)
        {
            sum=d_data[j-tid];

            for(int d=1; d<32; d*=2)
            {   
                T temp;
    
                if(std::is_same<T,FLOAT2>::value)
                    *(uint64_t*)&temp=__shfl_up_sync(-1,*(uint64_t*)&sum,d);
                else 
                    temp=__shfl_up_sync(-1,sum,d);
    
                if(tid>=d)
                    sum=sum+temp;
            }

            d_data[j-tid]=(float)(sum+pre_sum);

            if(std::is_same<T,FLOAT2>::value)
                *(uint64_t*)&sum=__shfl_sync(-1,*(uint64_t*)&sum,31);
            else 
                sum=__shfl_sync(-1,sum,31);
            pre_sum=pre_sum+sum;
        }

        if(j-tid>=0)
            sum=d_data[j-tid];

        for(int d=1; d<j; d*=2)
        {   
            T temp;

            if(std::is_same<T,FLOAT2>::value)
                *(uint64_t*)&temp=__shfl_up_sync(-1,*(uint64_t*)&sum,d);
            else 
                temp=__shfl_up_sync(-1,sum,d);

            if(tid>=d)
                sum=sum+temp;
      }

        if(j-tid>=0)
            d_data[j-tid]=(float)(sum+pre_sum);

    }//warpscan_backward
    
   template<typename T>
    __global__ void gpu_scansum_warp_v0(float *d_data, int length, int nbatch)
    {
        int tidx=threadIdx.x;
        int tidy=threadIdx.y;

        int n_per=blockDim.y*gridDim.x;

        for(int i=tidy+blockIdx.x*blockDim.y; i<nbatch; i+=n_per)
        {
            float *array=d_data+i*length;
            warpscan_forward<T>(tidx, array, length);
            warpscan_backward<T>(tidx,array, length);
        }
    }

性能优化

Grid和Threads运行配置对性能影响
nvprof 性能数据收集分析
Nsight Compute 性能数据收集分析
cudaEvent进行内核运行计时

环境#1运行时间结果

表中数字是每个kernel在GPU上的运行时间，单位毫秒
在这里插入图片描述

结论

warp版本整体效果最差，尤其是累加计算采用double和float2时候。一个原因是__shfl_up_sync调用时warp内有一定数量的线程等待不干活。warp版本对内存对齐不敏感，所有操作是在寄存器内进行，L2+L1 cache会一定程度上掩盖内存对齐问题!
gpu_smem_v1整体表现最好，但对数据内存对齐比较敏感，移除分支对性能提高有帮助
合适的grid和threads参数应该保证每个SM有一定数量的线程块运行，才能保证负载均衡，但也不要太大，会增加GPU调度时间
对shared memory的访问应该采用类似PADDING方式，尽可能避免bank conflict!
float2的版本明显比double版本性能高。GPU天生对单精度浮点支持最好

总结思考

使用shared memory,避免bank conflict对性能有较大影响，尤其对内存访问密集型操作
kernel运行布局配置(Grid & Block )对性能有影响，要保证调度给每个SM分配到足够多的线程块，有一定负载！
warp shuffles调用多次或是进行8Byte宽类型的移动操作，对性能有明显影响(可能不同体系结构结论会有区别)
对内存访问密集运算，输入数据在global memory的布局，会产生明显影响，应尽可能保证一定的对齐，对齐的内存访问会有明显性能提升
对7.0(Volta)以上计算能力GPU, warp同步不是显式的，需要调用__syncwarp
减少内层if语句
循环变量用int，避免用uint, 避免溢出问题！

weixin_42849849

关注

0
点赞
踩
0

收藏

觉得还不错? 一键收藏
0
评论
GPU实现Prefix Sum Scan (cumulative sum )算法

背景对给定NBATCH组数据float Array[NBATCH][LENGTH];对每组数据需要进行一次forward正向scan sum运算，然后进行backward反向一次scan sum。这种操作来源于地震道数据在偏移前的反假频滤波需要。其中NBATCH数量级可在10k左右，LENGTH数量级在10k左右。现有Thrust或CUB库没有提供适合多个数组的批量操作，对正反向结合也需要较复杂的逻辑表达。此外对于LENGTH较大的数组，需要考虑单精度数累加过程中精度丢失问题(rounding e
复制链接

扫一扫