Radix sort 基数排序

有关《GPU并行程序设计》(英文《 CUDA Programming A Developer’s Guide to Parallel Computing with GPUs》) 第六章 中基数排序，其中并行排序的多线程排序，由于没有具体较为详细的介绍，对于初次接触多线程的人略微困难。本文较为详细的介绍此多线程基数排序代码。

一：在解释多线程的代码之前，先解释下单线程的串行代码，如下：

__host__ void cpu_sort(u32 * const data,const u32 num_elements)
{
        static u32 cpu_tmp_0[NUM_ELEM];
        static u32 cpu_tmp_1[NUM_ELEM];

        for (u32 bit=0;bit<32;bit++)
        {
                u32 base_cnt_0 = 0;
                u32 base_cnt_1 = 0;

                for (u32 i=0; i<num_elements; i++)
                {
                        const u32 d = data[i];
                        const u32 bit_mask = (1 << bit);
                        if ( (d & bit_mask) > 0 )
                        {
                                cpu_tmp_1[base_cnt_1] = d;
                                base_cnt_1++;
                        }
                        else
                        {
                                cpu_tmp_0[base_cnt_0] = d;
                                base_cnt_0++;
                        }
                }
                // Copy data back to source - first the zero list
                 for (u32 i=0; i<base_cnt_0; i++)
                 {
                         data[i] = cpu_tmp_0[i];
                 }
                 // Copy data back to source - then the one list
                 for (u32 i=0; i<base_cnt_1; i++)
                 {
                         data[base_cnt_0+i] = cpu_tmp_1[i];
                 }
        }
}

串行代码执行过程，假设有 4 个整型数据元素的数组 data ，并创建两临时空间cpu_tmp_0，cpu_tmp_1 如下图：

只知考虑外循环bit=0时的循环情况，内层循环变量data中所有元素，执行过程如下图：

内存循环首先访问第一个元素4，判断4的末尾bit是0，所有将元素4存储到cpu_tmp_0的第一个位置；然后访问第二个元素3，判断3的末尾bit是1，所有将元素3存储到cpu_tmp_1的第一个位置；以此类推：以此方位第三、四个元素，分别存储到cpu_tmp_0、cpu_tmp_1的第二个位置。最后先将cpu_tmp_0的元素拷贝到data中对应位置中，在拷贝cpu_tmp_1中的元素。

当bit为其它值是，内层循环重复之前的计算过程。这个计算完成之后形成 一个有序序列！

二：如果理解的单线程的计算过程，那么多线程与单线程的 区别在于：

1、每个线程访问的元素不再连续的，而是步长为n，n为线程总数。

2、计算结束之后，形成了n个有序的序列。

先看下多线程代码，如下：

__device__ void radix_sort(u32 * const sort_tmp,const u32 num_lists,
			        const u32 num_elements,const u32 tid,
				u32 * const sort_tmp_0,u32 * const sort_tmp_1)
{
        // Sort into num_list, lists
        // Apply radix sort on 32 bits of data
        for (u32 bit=0;bit<32;bit++)
        {
                u32 base_cnt_0 = 0;
                u32 base_cnt_1 = 0;
                const u32 bit_mask = (1 << bit);
                for (u32 i=0; i<num_elements; i+=<strong>num_lists</strong>)
                {
                        const u32 elem = sort_tmp[i+tid];
                        if ( (elem & bit_mask) > 0 )
                        {
                                sort_tmp_1[base_cnt_1+tid] = elem;
                                base_cnt_1+=num_lists;
                        }
                        else
                        {
                                sort_tmp_0[base_cnt_0+tid] = elem;
                                base_cnt_0+=num_lists;
                        }
                }
        }
        // Copy data back to source - first the zero list
        for (u32 i=0; i<base_cnt_0; i+=num_lists)
        {
                sort_tmp[i+tid] = sort_tmp_0[i+tid];
        }
        // Copy data back to source - then the one list
        for (u32 i=0; i<base_cnt_1; i+=num_lists)
        {
                sort_tmp[base_cnt_0+i+tid] = sort_tmp_1[i+tid];
        }
       __syncthreads();
}

在多线程中变量 base_cnt_0，base_cnt_1为每个线程私有变量，每个线程的步长为 num_lists，num_lists在这里等于block中的线程数 。

在这里也只知考虑外循环bit=0时的循环情况，假设内层循环变量sort_tmp中共有如上所示16个元素，假设只有1个block，block中有4个线程，那么总共就是4个线程。那么 num_lists的值就是4，即：每个线程的步长是4.

如上图所示，4个线程每个线程处理一个元素，每次循环就会处理4个元素。共有16个元素，也就是循环4次。0号线程只处理索引为0、4、8、12的元素（红色块），同理，1号线程只处理索引是1、5、9、13（黄色快）；2号线程只处理索引为2、6、10、14的元素（蓝色快）；3号线程只处理索引为3、7、11、15的元素（绿色块）。

只看0号线程，红色块部分，处理过程，如下图：

0号线程只处理红色块的数据元素，且只能将sort_tmp中元素存放到，sort_tmp_0和sort_tmp_1中对应的红色位置内，最后在将sort_tmp_0和sort_tmp_1中的元素拷贝回sort_tmp。

同理，对于1、2、3号线程结果如下，

这是bit=0的结果，bit为其它值，处理过程也是一样的。最终将会得到4个有序的序列。即：红、黄、蓝、绿色块分别组成的序列。

至此基数排序结束。当然要想得到最后的一个完成的有序序列，则需后面章节讲到的-合并。

（完）