算法设计：基数排序
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 = ; bit < ; bit++)
{
u32 bit_mask = << bit, cnt0 = , cnt1 = ;
for (u32 i = ; i < num_elements; i ++)
{
u32 elem = sort_tmp[i];
if ((elem&bit_mask)>)
{
sort_tmp1[cnt1] = elem;
cnt1++;
}
else
{
sort_tmp[cnt0] = elem;
cnt0++;
}
}
for (u32 i = ; 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 = ; bit < ; bit++)
{
u32 bit_mask = << bit, cnt0 = , cnt1 = ;
for (u32 i = ; i < num_elements&&i + tid<num_elements; i += num_lists)
{
u32 elem = sort_tmp[i + tid];
if ((elem&bit_mask)>)
{
sort_tmp1[cnt1 + tid] = elem;
cnt1 += num_lists;
}
else
{
sort_tmp[cnt0 + tid] = elem;
cnt0 += num_lists;
}
}
for (u32 i = ; 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<<<,,>>>
即用内核函数的第三个参数指明动态数组大小，经过试验，在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[], *p4 = &p2[];
解释是，让一个指针指向共享内存的首地址，然后开动态空间，不过这奇葩的开法是错的，起码在CUDA 7.0里是不行的。
后来我又意识到，既然共享内存没用cudaMalloc开，而采用C方式，那么new会不会有用呢？我将代码换成：
extern __shared__ u32 sort_tmp[], sort_tmp1[];
int num1=,num2=
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] = ;//list_idx数组记录每个LIST的当前元素头
__syncthreads();
for (u32 i = ; 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 == ) { 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] = ;//list_idx数组记录每个LIST的当前元素头
__syncthreads();
for (u32 i = ; i < num_elements; i++)
{
u32 idx = list_idx[tid] * NUM_LISTS + tid;
u32 mid = NUM_LISTS >> ; //折半
if (idx<num_elements) elem = sort_tmp[idx];//各个线程取出元素，越界检查
else elem = inf;
reduction[tid] = elem; //构成临时归约数组
reduction_idx[tid] = tid;
__syncthreads(); //块内阻塞同步
while (mid != )
{
if (tid < mid) //屏蔽一半线程
{
u32 val2 = reduction[tid + mid];
if (reduction[tid] > val2) //对比两半线程
{
reduction[tid] = val2;
reduction_idx[tid] = tid + mid;
}
}
mid >>= ;//折半
__syncthreads();//块内阻塞同步,注意位置
}
if (tid == )
{
list_idx[reduction_idx[]]++; //元素头+1
data[i] = reduction[]; //写回显存
}
}
}