CUDA程序设计(三)

算法设计：基数排序

CUDA程序里应当尽量避免递归，因而在迭代排序算法里，基数排序通常作为首选。

1.1 串行算法实现

十进制位的基数排序需要考虑数位对齐问题，比较麻烦。通常实现的是二进制位的基数排序。

整体思路：与当前位做AND运算，按照0.....1的顺序重置序列，直到所有位迭代完毕。

sort_tmp数组作为基数桶，sort_tmp1作为辅助桶，存放当前位为1的数据。

__host__ void radix_sort(u32 *data,u32 *sort_tmp,u32 *sort_tmp1,u32 num_elements)
{
    for (u32 bit = 0; bit < 32; bit++)
    {
        u32 bit_mask = 1 << bit, cnt0 = 0, cnt1 = 0;
        for (u32 i = 0; i < num_elements; i ++)
        {
            u32 elem = sort_tmp[i];
            if ((elem&bit_mask)>0)
            {
                sort_tmp1[cnt1] = elem;
                cnt1++;
            }
            else
            {
                sort_tmp[cnt0] = elem;
                cnt0++;
            }
        }
        for (u32 i = 0; i < cnt1; i ++) sort_tmp[cnt0 + i] = sort_tmp1[i];
    }
}

 1.2 并行算法实现

基于数据分解的串改唯一注意点是，让相邻线程访问相邻数据，而不要让同一线程连续访问相邻数据。

经过多线程分解数据并行处理后，任何排序算法都会变成归并排序的中间状态。

__device__ void radix_sort(u32 *data,u32 *sort_tmp,u32 *sort_tmp1,u32 num_lists, u32 num_elements, u32 tid)
{
    for (u32 bit = 0; bit < 32; bit++)
    {
        u32 bit_mask = 1 << bit, cnt0 = 0, cnt1 = 0;
        for (u32 i = 0; i < num_elements&&i + tid<num_elements; i += num_lists)
        {
            u32 elem = sort_tmp[i + tid];
            if ((elem&bit_mask)>0)
            {
                sort_tmp1[cnt1 + tid] = elem;
                cnt1 += num_lists;
            }
            else
            {
                sort_tmp[cnt0 + tid] = elem;
                cnt0 += num_lists;
            }
        }
        for (u32 i = 0; i < cnt1; i+=num_lists) sort_tmp[cnt0 + i + tid] = sort_tmp1[i+tid];
    }
    __syncthreads();
}

1.3 性能分析

假设sort_tmp、sort_tmp1都是全局内存，且每个线程处理10个元素

那么R\W各：32*(10+5)=480次，每次500个T周期，这个时间是非常慢的。

共享内存与全局内存

2.1 共享内存机制

CUDA共享内存由线程块共享，默认连接着L1 Cache，因而访问有特别限制。

如果让一个线程连续访问相邻数据会怎么样？一个线程霸占着全部Cache，其它线程分不到Cache。

而这个线程后续数据还没有用到，却霸占着Cache的位置。其它线程分不到Cache，速度慢。

一旦__syncthreads后，需要等待最慢的线程结束，这样会导致Cache基本是废的。

这就是CUDA共享内存的 ”存储体冲突" 问题。无论是CPU还是GPU的Cache，都会出现这个问题。

罪魁祸首是基于数据分解的多线程算法模型。而CPU算法通常都是串行的，因而通常不是关注重点。

2.2 共享内存的使用方法

2.2.1 静态数组

开静态数组是基本手段，方法如下：

#define NUM_ELEMENTS XXXX
__shared__ u32 sort_tmp[NUM_ELEMENTS], sort_tmp1[NUM_ELEMENTS];

有趣的是，CUDA给__shared__设定的生存周期是整个线程块的周期，这意味着，

__shared__变量可以随地开，全局开也行，函数里开也行，不会转到栈空间去。

2.2.2 动态数组

CUDA早期的资料通常这样写着开动态数组的方法：

extern __shared__ u32 sort_tmp[], sort_tmp1[];
kernel_func<<<1,256,2333>>>

即用内核函数的第三个参数指明动态数组大小，经过试验，在CUDA 7.0中是无效的，目测官方已经废弃。

放弃的原因很简单，用统一的参数，只能开统一的大小，要是不同的大小呢？

大部分CUDA资料上几乎没有共享内存的指针申请法，唯一可追询的是这 http://blog.sina.com.cn/s/blog_5e8e35510100liz9.html

作者是这么做的：

extern __shared__ u32 sort_tmp[], sort_tmp1[];
u32 *p1 = sort_tmp, *p2 = sort_tmp1; 
u32 *p3 = &p1[2333], *p4 = &p2[2334];

解释是，让一个指针指向共享内存的首地址，然后开动态空间，不过这奇葩的开法是错的，起码在CUDA 7.0里是不行的。

后来我又意识到，既然共享内存没用cudaMalloc开，而采用C方式，那么new会不会有用呢？我将代码换成：

extern __shared__ u32 sort_tmp[], sort_tmp1[];
int num1=2333,num2=2334
u32 *p1 = sort_tmp, *p2 = sort_tmp1; 
p1 = new u32[num1], p2 = new u32[num2];

这回终于把动态共享内存开出来了。

2.3 全局内存机制

全局内存是CUDA最广泛存储体，由cudaMalloc申请，完全依附于显存，无权限进入Cache。

显存的访存周期长达500~600个T周期，为了没有Cache的缺陷，NVIDIA设计了线程束访存机制。

与共享内存的数据排布类似，该机制让相邻线程访问相邻数据，最小限制单位是half-warpSize(16个线程)

只要相邻的16个线程访问相邻的全局内存，就可以获得最大128字节的一步预读。

归并

3.1 并行合并

一共进行N轮推选，每轮中，各个线程返回元素序列头，决出最值。

__device__ void merge_parallel(u32 *data,u32 *sort_tmp,u32 num_elements,u32 tid)
{
    __shared__ u32 min_value, min_tid;
    __shared__ u32 list_idx[NUM_LISTS];//共享内存，访问越频繁，Cache利用率越高
    u32 elem;list_idx[tid] = 0;//list_idx数组记录每个LIST的当前元素头
    __syncthreads();
    for (u32 i = 0; i < num_elements; i++)
    {
        u32 idx = list_idx[tid] * NUM_LISTS + tid;
        //注意：共享内存的存放方式
        //线程的下一个元素需要跳跃NUM_LISTS单位，tid则决定着是哪个LIST
        if (idx<num_elements) elem = sort_tmp[idx];//各个线程取出元素，越界检查
        else elem = inf;
        if (tid == 0) { min_value = min_tid = inf; } //初始化
        __syncthreads(); //块内阻塞同步
        atomicMin(&min_value, elem); //块内原子求最小值
        __syncthreads();//块内阻塞同步
        //线程检查：如果块内最小值是自己提供的，则上报
        //二次检查：如果有多个上报邀功的，则取最小tid的
        if (min_value == elem) atomicMin(&min_tid, tid);
        __syncthreads();//块内阻塞同步：防止未决出最小值，就向下执行
        if (min_tid == tid)
        {
            list_idx[tid]++; //元素头+1
            data[i] = elem; //写回显存
        }
    }
}

3.2 并行二分归约

CUDA 1.2以下版本没有原子函数，所以得使用另一种既能找出最值，又能避免访存冲突的方法。

N个数求最值可以通过并行在$log(N)$时间内完成。

每轮中，将数据一分为二，前一半与后一半比较，将最值写回前一半。完成二分。

当然$log(N)$并不会载入史册，因为这是个错误的复杂度分析。

尽管仅需要$log(N)$轮，但每轮中，比较的分组是不可能完全并行的。

CUDA中理论最大并行线程是2048。如果有10000个数据，那么第一轮需要5000组比较：

CPU串行执行要循环5000次。

GPU并行也要循环：[5000/2048]=2次

令人惊讶的是，CPU串行归约第一轮就需要循环5000次，还不如不归约了。

二分归约时间直接依赖于同时并行量，并行量越大，效率越高。

反之，并行量越小，效率越低。在串行中，甚至退化成了负效率。

__device__ void merge_parallel2(u32 *data, u32 *sort_tmp, u32 num_elements, u32 tid)
{
    __shared__ u32 list_idx[NUM_LISTS];//共享内存，访问越频繁，Cache利用率越高
    __shared__ u32 reduction[NUM_LISTS], reduction_idx[NUM_LISTS];
    u32 elem; list_idx[tid] = 0;//list_idx数组记录每个LIST的当前元素头
    __syncthreads();
    for (u32 i = 0; i < num_elements; i++)
    {
        u32 idx = list_idx[tid] * NUM_LISTS + tid;
        u32 mid = NUM_LISTS >> 1; //折半
        if (idx<num_elements) elem = sort_tmp[idx];//各个线程取出元素，越界检查
        else elem = inf;
        reduction[tid] = elem; //构成临时归约数组
        reduction_idx[tid] = tid;
        __syncthreads(); //块内阻塞同步
        while (mid != 0)
        {
            if (tid < mid) //屏蔽一半线程
            {
                u32 val2 = reduction[tid + mid];
                if (reduction[tid] > val2)  //对比两半线程
                {
                    reduction[tid] = val2;
                    reduction_idx[tid] = tid + mid;
                }
            }
            mid >>= 1;//折半
            __syncthreads();//块内阻塞同步,注意位置
        }
        if (tid == 0)
        {
            list_idx[reduction_idx[0]]++; //元素头+1
            data[i] = reduction[0]; //写回显存
        }
    }
}