作者 | Derrick 编辑 | 自动驾驶之心
原文链接：https://zhuanlan.zhihu.com/p/720276184
点击下方卡片，关注“自动驾驶之心”公众号
戳我-> 领取自动驾驶近15个方向学习路线
本文只做学术分享，如有侵权，联系删文
一、前言
排序（sort），是计算机内经常进行的一种操作，其目的是将一组“无序”的记录序列调整为“有序”的记录序列，是各领域中广泛应用的核心算法之一。受限于算法本身的限制（所有基于比较的排序算法复杂度下限为O(n)，对于一组超大量数据单纯的cpu排序所需时间过长，通常1000万的数据排序就需要1秒左右，更大量的数据排序则需要分钟级的运行时间。因此，可以通过并行排序操作来进一步加速排序算法。
二、基于cpu的排序算法
经过数十年的发展，目前的排序算法已有基数排序、冒泡排序、选择排序、插入排序、希尔排序、堆排序、归并排序、快速排序等。显然，我们希望在已有的高效排序算法上进行进一步的优化，因些我们忽略冒泡排序、选择排序、希尔排序这一类指数级复杂度的算法，也忽略基数排序这一类不适于于浮点型变量的算法，而在堆排序、归并排序、快速排序这些O(n)算法中，直观的我们能发现归并排序、快速排序都是基于二分的策略，所以可以自然地同时对左右两半同时进行排序；而归并的过程也可以通过直接查找每个元素在合并数组中的位置来契合并行化，但快速排序中划分的过程无法并行实现（至少我没有想到，欢迎各路大神提出好的思路）；因此归并排序是最适于并行化的排序算法。
然而，这些已有的高效排序算法通常只在数据量较大时才有着更高的效率，在小数据时反而常数更小的选择排序、插入排序等有着更出色的性能；如中的std::sort中在长度小于16或32时不再继续划分，而是调用插入排序来实现；因此在基于cuda的排序中也需要在小数据上选择另一种算法以达到更出色的性能。在此，我们选择并不常见、但复杂度O(nn)相对较低的双调排序来实现，自然是因为它能很好的并行化。
因此我们选择双调排序和归并排序来实现CUDA上的排序操作，如需详细了解这两个算法请查阅相关资料。
三、基于cuda的排序算法的实现
1. 双调排序的并行化
我们暂定一个线程块的线程数为 256（tile_size），因为双调排序中每次需要进行一次比较来判断是否需要交换，因此每一个线程块处理连续的tile_size x 2个数据（以下称为一个数据段），每个线程负责一次比较操作。因此该核函数线程块大小为tile_size，网格大小为总共的数据段数量。
#define tile_size 256u
#define log_tile_size 8u
unsigned int b = n + (tile_size << 1) - 1 >> log_tile_size + 1, s = tile_size << 1;
binoticsort_gpu<<<b, s >> 1>>>(data, n);首先先将需要的数据载入共享内存中，由于每一个线程块处理tile_size x 2个数据，所以一个线程块负责的数据段起始为“线程块索引 x tile_size x 2”；每一个线程将相隔tile_size的两个数据载入共享内存中，超过数据规模的部分直接填充浮点数的最大值，这样一个线程束访问一段连续的地址以达到更高的访存效率；还需添加一个线程块内同步函数，确保所有数据都已载入共享内存后再开始排序。
__shared__ double buffer[tile_size << 1];
unsigned int d = blockIdx.x << log_tile_size + 1, i = threadIdx.x;
data += d; n -= d;
buffer[i] = i < n ? data[i] : DBL_MAX;
buffer[i + blockDim.x] = i + blockDim.x < n ? data[i + blockDim.x] : DBL_MAX;
__syncthreads();然后需要将随机序列转换为双调序列，在这个过程中需要交换操作，因此需要先定义一个设备端的交换函数，该交换函数可定义为模板函数以支持不同类型的数据交换。
// auxiliary functions
template<typename type>__device__ void swap(type& v1, type& v2)
{
type t = v1;
v1 = v2;
v2 = t;
}
unsigned int s, t;
for (unsigned int k = 2; k <= blockDim.x; k <<= 1)
for (unsigned int p = k >> 1; p; p >>= 1)
{
s = (i & -p) << 1 | i & p - 1;
t = s | p;
if (s & k ? buffer[s] < buffer[t] : buffer[s] > buffer[t])
swap(buffer[s], buffer[t]);
__syncthreads();
}
该过程转化为代码的操作与上述合并过程中的内循环完全一致。
for (unsigned int p = blockDim.x; p; p >>= 1)
{
s = (i & -p) << 1 | i & p - 1;
t = s | p;
if (buffer[s] > buffer[t])
swap(buffer[s], buffer[t]);
__syncthreads();
}最后将该数据段的有效排序结果写回全局内存。
if (i < n)
data[i] = buffer[i];
if (i + blockDim.x < n)
data[i + blockDim.x] = buffer[i + blockDim.x];最终分块双调排序的代码如下：
#include <float.h>
// binotic sort
__global__ void binoticsort_gpu(double* data, unsigned int n)
{
__shared__ double buffer[tile_size << 1];
unsigned int d = blockIdx.x << log_tile_size + 1, i = threadIdx.x;
data += d; n -= d;
buffer[i] = i < n ? data[i] : DBL_MAX;
buffer[i + blockDim.x] = i + blockDim.x < n ? data[i + blockDim.x] : DBL_MAX;
__syncthreads();
unsigned int s, t;
for (unsigned int k = 2; k <= blockDim.x; k <<= 1)
for (unsigned int p = k >> 1; p; p >>= 1)
{
s = (i & -p) << 1 | i & p - 1;
t = s | p;
if (s & k ? buffer[s] > buffer[t] : buffer[s] < buffer[t])
swap(buffer[s], buffer[t]);
__syncthreads();
}
for (unsigned int p = blockDim.x; p; p >>= 1)
{
s = (i & -p) << 1 | i & p - 1;
t = s | p;
if (buffer[s] > buffer[t])
swap(buffer[s], buffer[t]);
__syncthreads();
}
if (i < n)
data[i] = buffer[i];
if (i + blockDim.x < n)
data[i + blockDim.x] = buffer[i + blockDim.x];
}经过该分块双调排序后，规模为n的数据被划分为多个长为 tile_size x 2的数据段，每个数据段内的 tile_size x 2个数据已经被成功排序，下面只需要逐步将这些有序数据段合并即可。
2. 归并操作的并行化
有序数据段合并操作即像归并排序中的归并一样，但是需要对其做一定的修改以更好地适用于并行化。通常的合并操作需从头逐个比较两个数组中的元素，取较小者置于结果数组中——但是这个循环是具有数据依赖性的，故这种相互之间不独立的操作无法进行并行化处理。因此我们需要用另一种方法来进行合并操作，我们直接计算两个待合并数组中每个元素在结果数组中应当存在的位置，然后把该元素置于该位置即可。
显然对于其中一待合并数组中比某元素小的元素个数必然是该元素的索引，因此只需找到另一待合并数组中小于（或小于等于）该元素的元素个数即可，注意这儿必须要一个数组用小于而另一个数组用小于等于，否则两个相同的数字就会落在同一位置。在此直接引用中的std::lower_bound和std::upper_bound中的算法来实现，只是需要将其编译为设备端函数。
// auxiliary functions
template<typename type>__device__ unsigned int lower_bound(const type* data, unsigned int n, const type& value)
{
unsigned int len = n, half;
const type* start = data, * mid;
while (len)
{
half = len >> 1;
mid = data + half;
if (*mid < value)
{
data = mid + 1;
len = len - half - 1;
}
else
len = half;
}
return data - start;
}
template<typename type>__device__ unsigned int upper_bound(const type* data, unsigned int n, const type& value)
{
unsigned int len = n, half;
const type* start = data, *mid;
while (len)
{
half = len >> 1;
mid = data + half;
if (value < *mid)
len = half;
else
{
data = mid + 1;
len = len - half - 1;
}
}
return data - start;
}
mergedirect_gpu<<<dim3(s >> log_tile_size, std::min(b, 32768u), b + 32767 >> 15), tile_size>>>(data, n, s, tmp);
// merge directly
__global__ void mergedirect_gpu(const double* data, unsigned int n, unsigned int s, double* result)
{
unsigned int d1 = (blockIdx.z * gridDim.y + blockIdx.y) * s << 1, d2 = d1 + s, i = d1 + blockIdx.x * blockDim.x + threadIdx.x;
if (i < n)
result[d2 < n ? i + lower_bound(data + d2, min(s, n - d2), data[i]) : i] = data[i];
if (s + i < n)
result[i + upper_bound(data + d1, s, data[s + i])] = data[s + i];
}然后将这两个函数整合成一个排序函数即可，该函数接受两个参数：指向数据首地址的指针和数据规模，返回指向已排序数据首地址的指针，其中的指针均指向设备端。该函数首先需要申请一块与输入数据等长的临时数组用于存放合并的结果，然后先进行一次双调排序将输入数据分块排序为多个长为 tile_size x 2的有序数据段，再来一个循环将这些有序数据段逐轮合并直至数据段长度大于数据规模，最后释放临时数组空间并返回结果。
double* sortdirect_gpu(double* data, unsigned int n)
{
double* tmp;
cudaMalloc(&tmp, sizeof(double) * n);
unsigned int b = n + (tile_size << 1) - 1 >> log_tile_size + 1, s = tile_size << 1;
binoticsort_gpu<<<b, s >> 1>>>(data, n);
for (b = b + 1 >> 1; s < n; s <<= 1, b = b + 1 >> 1)
{
mergedirect_gpu<<<dim3(s >> log_tile_size, std::min(b, 32768u), b + 32767 >> 15), tile_size>>>(data, n, s, tmp);
std::swap(data, tmp);
}
cudaDeviceSynchronize();
cudaFree(tmp);
return data;
}至此，我们已初步实现了基于CUDA的排序算法，我们选择了不同的数据规模在RTX A6000上对它的性能进行测试，测试代码如下：
// test.cu
#include <cuda.h>
#include <cuda_runtime.h>
#include <device_launch_parameters.h>
#include <iostream>
#include <chrono>
#include <algorithm>
int main(int argc, char** argv)
{
cudaSetDevice(0);
cudaFree(0);
unsigned int test[] = { 1025u, 65537u, 1048577u, 16777217u, 134217729u, 1073741825u };
const unsigned int runtime = 3;
std::chrono::high_resolution_clock::time_point start, end;
unsigned long long elasped[3] = { 0ull };
for (unsigned int t = 0; t < sizeof(test) / sizeof(unsigned int); ++t)
{
unsigned int n = test[t];
double* original = new double[n];
double* data = new double[n];
elasped[0] = elasped[1] = elasped[2] = 0ull;
srand(time(nullptr));
for (int i = 0; i < runtime; ++i)
{
for (unsigned int i = 0; i < n; ++i)
original[i] = data[i] = (rand() * 2.0 / RAND_MAX - 1.0) * rand();
// cpu sort 1
start = std::chrono::high_resolution_clock::now();
qsort(data, n, sizeof(double), [](const void* left, const void* right)->int {double tmp = *(double*)left - *(double*)right; if (tmp > 0) return 1; else if (tmp < 0) return -1; else return 0; });
end = std::chrono::high_resolution_clock::now();
elasped[0] += std::chrono::duration_cast<std::chrono::nanoseconds>(end - start).count();
// cpu sort 2
memcpy(data, original, sizeof(double) * n);
start = std::chrono::high_resolution_clock::now();
std::sort(data, data + n);
end = std::chrono::high_resolution_clock::now();
elasped[1] += std::chrono::duration_cast<std::chrono::nanoseconds>(end - start).count();
// gpu sort
double* pd_data;
cudaMalloc(&pd_data, sizeof(double) * n);
cudaMemcpy(pd_data, original, sizeof(double) * n, cudaMemcpyHostToDevice);
start = std::chrono::high_resolution_clock::now();
pd_data = sortdirect_gpu(pd_data, n);
end = std::chrono::high_resolution_clock::now();
elasped[2] += std::chrono::duration_cast<std::chrono::nanoseconds>(end - start).count();
cudaMemcpy(data, pd_data, sizeof(double) * n, cudaMemcpyDeviceToHost);
cudaFree(pd_data);
}
printf("data number: %u\n", n);
printf("cpu merge sort cost %.5f ms\n", elasped[0] / 1e6 / runtime);
printf("cpu quick sort cost %.5f ms\n", elasped[1] / 1e6 / runtime);
printf("gpu sort cost %.5f ms\n", elasped[2] / 1e6 / runtime);
printf("-------------------------------------------\n");
delete[] data;
delete[] original;
}
return 0;
}我们先不加任何优化进行编译
nvcc -o ./benchmark -g ./test.cu && ./benchmark其运行结果如下（只统计排序所消耗的时间，忽略数据传输的耗时）：
然后我们开启O2优化进行编译
nvcc -O2 -o ./benchmark_O2 -g ./test.cu && ./benchmark_O2其运行结果如下（只统计排序所消耗的时间，忽略数据传输的耗时）：
我们再加上数据传输耗时，其耗时如下：
可以发现在不计数据传输开销时，从数万到10亿的数据规模下，该基于CUDA的排序算法相较于单线程cpu排序算法而言都能获得数十乃至百倍的加速比；若算上数据传输的耗时，数据规模较小时其性能不如cpu排序算法，但在数据规模达100万及以上仍有着20倍左右的加速比（相较于更快的std::sort）。
四. 基于cuda的排序算法的改进
忽略数据传输的耗时，我们已经实现了约60倍的加速比（相较于更快的std::sort），但是仍能注意到该排序算法中合并操作是主要的耗时步骤，而其主要的耗时步骤则在查找每个元素在另一数据段中的位置，因为当数据段长度较大时，在很大的地址区间内进行二分查找是一个很费时的操作，这既不利于缓存（cache）的命中又不利于线程束（warp）的合并访问（coalesced access）。因此我们试图对需要合并的数据段进行进一步的分块，然后分别将两个数据段中对应的块进行合并，以进一步提升该排序算法的性能。
1. 分块合并
首先我们需要考虑如何对数据段进行合理的分块，因为必须要保证待合并的两个数据段第一块中的所有元素均不大于第二块中的所有元素，这样才能在第一对块合并前确认第二对块合并后的位置以保证并行性。最直观的方法为将第一个数据段直接分块，然后在第二个数据段中找到对应的块，可通过在第二个数据段中查找第一个数据段各块首元素的位置，按这些位置便可将第二个数据段分块。
这样的分块策略对于随机数据来说能取得不错的效果，但如果两个数据段中的元素分布不均匀，甚至第一个数据段中的所有元素都小于（或大于）第二个数据段，该分块策略就不能进行有效的切分，也不能达到良好的性能。因此除了第一个数据段各块首元素外，还用第二个数据段等分后每块首元素来进行切分，且同时对两个数据段进行切分以得到对应的块。
接下来仍然是最令人头疼（甚至是头秃）的时刻，就是把上述流程转换为核函数的代码。最后一步必然是将对应块进行合并，而对于切分操作有两种策略：一是先取出所有块首元素，排序后再用来切分两个数据段；二是先用所有块首元素来切分两个数据段，再将切分位点进行排序。显然前者需要3步而后者只需2步，鉴于核函数启动的耗时我们选择后者。
第一步是用所有块首元素来切分两个数据段，为此需要有一个rank数组来保存每个首元素的切分位点。我们让一个线程负责计算一个块首元素在两个数据段中的切分点，在当前数据段中的切分位点可直接由块索引乘以块大小得到，在另一数据段中的切分位点由二分查找得到。rank数组中依次保存第一、二个数据段的块首元素在第一、二个数据段中的切分位点，
__global__ void rank_gpu(const double* data, unsigned int n, unsigned int s, unsigned int* rank)
{
unsigned int pl = s >> log_tile_size, t = blockIdx.x * blockDim.x + threadIdx.x, k = t / (pl << 1), i = t & (pl << 1) - 1;
data += k * s << 1; n -= k * s << 1; rank += k * pl << 2;
unsigned int pr = min(n - s, s) + tile_size - 1 >> log_tile_size, p = pl + pr;
if (n <= s || i >= p)
return;
rank[i] = i & pl ? lower_bound(data, s, data[i << log_tile_size]) : i << log_tile_size;
rank[p + i] = i & pl ? (i ^ pl) << log_tile_size : lower_bound(data + s, min(n - s, s), data[i << log_tile_size]);
}第二步将所有数据段的两半切分位点进行排序，实际也是两半有序数组的合并，故其原理与直接合并操作相同，由于每个数据段可能会有超过 tile_size 个块，所以需要一个线程块组来负责一次合并；这儿我们选择让一个线程块组来负责一对数据段中切分位点的合并（两次合并），因为切分位点的总个数与这对数据段的长度相关（最后一对数据段中第二个数据段可能是不满的）。所有线程块组仍然排列为二维，因为网格大小第二、三维上限是-1。
首先仍然是数据段包含的块数量、、p及 data、rank、split数组偏移量的计算，除数据段组别k和当前线程负责的数据索引i分别由线程块索引的第二、三维和线程块索引的第一维、线程索引得到外，其余计算与第一步完全一致，若偏移量大于总数据规模则直接返回。最后在判断数据索引是否有效（i< 块数量）后，利用二分查找将切分位点进行排序，结果写入split数组中。
__global__ void ranksort_gpu(const unsigned int* rank, unsigned int n, unsigned int s, unsigned int* split)
{
unsigned int pl = s >> log_tile_size, k = blockIdx.z * gridDim.y + blockIdx.y, i = blockIdx.x * blockDim.x + threadIdx.x;
if (k * s << 1 >= n)
return;
n -= k * s << 1; rank += k * pl << 2; split += k * pl << 2;
unsigned int pr = min(n - s, s) + tile_size - 1 >> log_tile_size, p = pl + pr;
if (i < pl)
{
split[i + lower_bound(rank + pl, pr, rank[i])] = rank[i];
split[p + i + lower_bound(rank + p + pl, pr, rank[p + i])] = rank[p + i];
}
if (i < pr)
{
split[i + upper_bound(rank, pl, rank[pl + i])] = rank[pl + i];
split[p + i + upper_bound(rank + p, pl, rank[p + pl + i])] = rank[p + pl + i];
}
}第三步则是根据每组数据段的切分位点进行分块合并，每个线程块负责一对对应块的合并，因此每对数据段的合并需要一组线程块来完成。
__global__ void mergesegment_gpu(const double* data, const unsigned int* split, unsigned int n, unsigned int s, double* result)
{
__shared__ double databuffer[tile_size << 1], resultbuffer[tile_size << 1];
unsigned int pl = s >> log_tile_size, k = blockIdx.z * gridDim.y + blockIdx.y;
if (k * s << 1 >= n)
return;
data += k * s << 1; result += k * s << 1; n -= k * s << 1; split += k * pl << 2;
unsigned int pr = min(n - s, s) + tile_size - 1 >> log_tile_size, p = pl + pr;
if (n <= s)
{
unsigned int i = blockIdx.x * blockDim.x + threadIdx.x;
if (i < n)
result[i] = data[i];
}
else if (blockIdx.x < p)
{
unsigned int fstart = split[blockIdx.x], lstart = split[p + blockIdx.x];
unsigned int fend = blockIdx.x + 1 < p ? split[blockIdx.x + 1] : s, lend = blockIdx.x + 1 < p ? split[p + blockIdx.x + 1] : min(n - s, s);
unsigned int flen = fend - fstart, llen = lend - lstart;
if (threadIdx.x < flen)
databuffer[threadIdx.x] = data[fstart + threadIdx.x];
if (threadIdx.x < llen)
databuffer[flen + threadIdx.x] = data[s + lstart + threadIdx.x];
__syncthreads();
if (threadIdx.x < flen)
resultbuffer[threadIdx.x + lower_bound(databuffer + flen, llen, databuffer[threadIdx.x])] = databuffer[threadIdx.x];
if (threadIdx.x < llen)
resultbuffer[threadIdx.x + upper_bound(databuffer, flen, databuffer[flen + threadIdx.x])] = databuffer[flen + threadIdx.x];
__syncthreads();
if (threadIdx.x < flen)
result[fstart + lstart + threadIdx.x] = resultbuffer[threadIdx.x];
if (threadIdx.x < llen)
result[fstart + lstart + flen + threadIdx.x] = resultbuffer[flen + threadIdx.x];
}
}最后将上面三个函数及之前的分块双调排序函数组合成最终的排序函数，这一改进版的排序函数中除了用于存放临时合并结果的数组外，还额外需要两个索引数组用于存放切分位点和排序后的切分位点，每个索引数组的长度应是所有数据的总分块数 x 2，因为每个数据段还要被同组另一数据段中的块首元素切分。仍先进行一次双调排序进行分块排序，然后逐轮分块合并，最后释放临时数组空间并返回结果。注意分块合并中第一步的总线程数是固定的，即为总分块数，因为一个线程负责索引数组中两个索引的计算；第二步中一组线程块负责一对数据段中切分位点的排序，因此线程块组的总数即为长为2s的数据段对数（即为网格的第二、三维大小），而每个线程负责四个索引的排序，所以一对数据段中切分位点的排序总共需要的线程数恰等于每个数据段的初始分块数（即为网格的第一维大小和线程块大小）；第三步中一组线程块负责一对数据段的分块合并，所以其线程块组的总数与第二步相同，但若存在额外的不足一个数据段的剩余数据（所有数据的总分块数是偶数时）则需要一个额外的线程块组来进行额外数据的拷贝，而每组线程块需要的线程块数为每个数据段的分块数，恰等于其初始分块数 x 2（即为网格的第一维大小）。
double* sortsegment_gpu(double* data, unsigned int n)
{
double* tmp;
cudaMalloc(&tmp, sizeof(double) * n);
unsigned int* rank, index_len = (n + tile_size - 1 >> log_tile_size) << 1;
cudaMalloc(&rank, sizeof(unsigned int) * index_len << 1);
unsigned int* split = rank + index_len;
unsigned int b = n + (tile_size << 1) - 1 >> log_tile_size + 1, s = tile_size << 1, m = n + tile_size - 1 >> log_tile_size;
binoticsort_gpu<<<b, s >> 1>>>(data, n);
for (b = b + 1 >> 1; s < n; s <<= 1, b = b + 1 >> 1)
{
unsigned int p = s + tile_size - 1 >> log_tile_size, k = (n - 1) / s + 1 >> 1;
rank_gpu<<<m + tile_size - 1 >> log_tile_size, tile_size>>>(data, n, s, rank);
ranksort_gpu<<<dim3(p + tile_size - 1 >> log_tile_size, std::min(k, 32768u), k + 32767 >> 15), tile_size>>>(rank, n, s, split);
if (!((n - 1) / s & 1))
++k;
mergesegment_gpu<<<dim3(p << 1, std::min(k, 32768u), k + 32767 >> 15), tile_size>>>(data, split, n, s, tmp);
std::swap(data, tmp);
}
cudaDeviceSynchronize();
cudaFree(tmp);
cudaFree(rank);
return data;
}至此我们已经实现了基于分块合并的改进版 CUDAsort，下面我们对其性能进行评估。
2. 性能测试
我们用相同的数据规模对该改进版进行性能测试，测试结果如下：
我们发现数据规模在百万及以下时该改进版的性能反而有所下降，这是因为在数据规模较小时GPU的缓存能较好地缓解对全局内存频繁访问的时延，而此时分块合并法由于有更多的核函数启动会消耗更多的时间，故其性能反而不如直接合并法；但当数据规模在千万及以上时该改进版展现了更优异的性能，10亿双精度浮点型数据的排序仅需约1.1S，相较于初始版本有40%的性能提升。
五. 总结
我们基于双调排序和归并排序实现了GPU上的排序算法，若忽略数据传输的耗时，改进版 CUDAsort在数据规模较大时能达到百倍左右的加速比，能在1s左右完成10亿数据的排序。但是作者发现在tile_size设为512或1024时，该排序算法似乎会出错，思索半天未果，在此欢迎各路大神发表看法，包括但不限于出错的原因、更高性能的排序算法等，本人由于并不是计算机专业的，水平有限，因此只能想到这一种GPU上的高效排序算法了。
最后贴上完整的代码：
#include <cuda.h>
#include <cuda_runtime.h>
#include <device_launch_parameters.h>
#include <float.h>
#define tile_size 256u
#define log_tile_size 8u
// auxiliary functions
template<typename type>__device__ void swap(type& v1, type& v2)
{
type t = v1;
v1 = v2;
v2 = t;
}
template<typename type>__device__ unsigned int lower_bound(const type* data, unsigned int n, const type& value)
{
unsigned int len = n, half;
const type* start = data, * mid;
while (len)
{
half = len >> 1;
mid = data + half;
if (*mid < value)
{
data = mid + 1;
len = len - half - 1;
}
else
len = half;
}
return data - start;
}
template<typename type>__device__ unsigned int upper_bound(const type* data, unsigned int n, const type& value)
{
unsigned int len = n, half;
const type* start = data, *mid;
while (len)
{
half = len >> 1;
mid = data + half;
if (value < *mid)
len = half;
else
{
data = mid + 1;
len = len - half - 1;
}
}
return data - start;
}
// binotic sort
__global__ void binoticsort_gpu(double* data, unsigned int n)
{
__shared__ double buffer[tile_size << 1];
unsigned int d = blockIdx.x << log_tile_size + 1, i = threadIdx.x;
data += d; n -= d;
buffer[i] = i < n ? data[i] : DBL_MAX;
buffer[i + blockDim.x] = i + blockDim.x < n ? data[i + blockDim.x] : DBL_MAX;
__syncthreads();
unsigned int s, t;
for (unsigned int k = 2; k <= blockDim.x; k <<= 1)
for (unsigned int p = k >> 1; p; p >>= 1)
{
s = (i & -p) << 1 | i & p - 1;
t = s | p;
if (s & k ? buffer[s] < buffer[t] : buffer[s] > buffer[t])
swap(buffer[s], buffer[t]);
__syncthreads();
}
for (unsigned int p = blockDim.x; p; p >>= 1)
{
s = (i & -p) << 1 | i & p - 1;
t = s | p;
if (buffer[s] > buffer[t])
swap(buffer[s], buffer[t]);
__syncthreads();
}
if (i < n)
data[i] = buffer[i];
if (i + blockDim.x < n)
data[i + blockDim.x] = buffer[i + blockDim.x];
}
// merge directly
__global__ void mergedirect_gpu(const double* data, unsigned int n, unsigned int s, double* result)
{
unsigned int d1 = (blockIdx.z * gridDim.y + blockIdx.y) * s << 1, d2 = d1 + s, i = d1 + blockIdx.x * blockDim.x + threadIdx.x;
if (i < n)
result[d2 < n ? i + lower_bound(data + d2, min(s, n - d2), data[i]) : i] = data[i];
if (s + i < n)
result[i + upper_bound(data + d1, s, data[s + i])] = data[s + i];
}
double* sortdirect_gpu(double* data, unsigned int n)
{
double* tmp;
cudaMalloc(&tmp, sizeof(double) * n);
unsigned int b = n + (tile_size << 1) - 1 >> log_tile_size + 1, s = tile_size << 1;
binoticsort_gpu<<<b, s >> 1>>>(data, n);
for (b = b + 1 >> 1; s < n; s <<= 1, b = b + 1 >> 1)
{
mergedirect_gpu<<<dim3(s >> log_tile_size, std::min(b, 32768u), b + 32767 >> 15), tile_size>>>(data, n, s, tmp);
std::swap(data, tmp);
}
cudaDeviceSynchronize();
cudaFree(tmp);
return data;
}
// merge in segments
__global__ void rank_gpu(const double* data, unsigned int n, unsigned int s, unsigned int* rank)
{
unsigned int pl = s >> log_tile_size, t = blockIdx.x * blockDim.x + threadIdx.x, k = t / (pl << 1), i = t & (pl << 1) - 1;
data += k * s << 1; n -= k * s << 1; rank += k * pl << 2;
unsigned int pr = min(n - s, s) + tile_size - 1 >> log_tile_size, p = pl + pr;
if (n <= s || i >= p)
return;
rank[i] = i & pl ? lower_bound(data, s, data[i << log_tile_size]) : i << log_tile_size;
rank[p + i] = i & pl ? (i ^ pl) << log_tile_size : lower_bound(data + s, min(n - s, s), data[i << log_tile_size]);
}
__global__ void ranksort_gpu(const unsigned int* rank, unsigned int n, unsigned int s, unsigned int* split)
{
unsigned int pl = s >> log_tile_size, k = blockIdx.z * gridDim.y + blockIdx.y, i = blockIdx.x * blockDim.x + threadIdx.x;
if (k * s << 1 >= n)
return;
n -= k * s << 1; rank += k * pl << 2; split += k * pl << 2;
unsigned int pr = min(n - s, s) + tile_size - 1 >> log_tile_size, p = pl + pr;
if (i < pl)
{
split[i + lower_bound(rank + pl, pr, rank[i])] = rank[i];
split[p + i + lower_bound(rank + p + pl, pr, rank[p + i])] = rank[p + i];
}
if (i < pr)
{
split[i + upper_bound(rank, pl, rank[pl + i])] = rank[pl + i];
split[p + i + upper_bound(rank + p, pl, rank[p + pl + i])] = rank[p + pl + i];
}
}
__global__ void mergesegment_gpu(const double* data, const unsigned int* split, unsigned int n, unsigned int s, double* result)
{
__shared__ double databuffer[tile_size << 1], resultbuffer[tile_size << 1];
unsigned int pl = s >> log_tile_size, k = blockIdx.z * gridDim.y + blockIdx.y;
if (k * s << 1 >= n)
return;
data += k * s << 1; result += k * s << 1; n -= k * s << 1; split += k * pl << 2;
unsigned int pr = min(n - s, s) + tile_size - 1 >> log_tile_size, p = pl + pr;
if (n <= s)
{
unsigned int i = blockIdx.x * blockDim.x + threadIdx.x;
if (i < n)
result[i] = data[i];
}
else if (blockIdx.x < p)
{
unsigned int fstart = split[blockIdx.x], lstart = split[p + blockIdx.x];
unsigned int fend = blockIdx.x + 1 < p ? split[blockIdx.x + 1] : s, lend = blockIdx.x + 1 < p ? split[p + blockIdx.x + 1] : min(n - s, s);
unsigned int flen = fend - fstart, llen = lend - lstart;
if (threadIdx.x < flen)
databuffer[threadIdx.x] = data[fstart + threadIdx.x];
if (threadIdx.x < llen)
databuffer[flen + threadIdx.x] = data[s + lstart + threadIdx.x];
__syncthreads();
if (threadIdx.x < flen)
resultbuffer[threadIdx.x + lower_bound(databuffer + flen, llen, databuffer[threadIdx.x])] = databuffer[threadIdx.x];
if (threadIdx.x < llen)
resultbuffer[threadIdx.x + upper_bound(databuffer, flen, databuffer[flen + threadIdx.x])] = databuffer[flen + threadIdx.x];
__syncthreads();
if (threadIdx.x < flen)
result[fstart + lstart + threadIdx.x] = resultbuffer[threadIdx.x];
if (threadIdx.x < llen)
result[fstart + lstart + flen + threadIdx.x] = resultbuffer[flen + threadIdx.x];
}
}
double* sortsegment_gpu(double* data, unsigned int n)
{
double* tmp;
cudaMalloc(&tmp, sizeof(double) * n);
unsigned int* rank, index_len = (n + tile_size - 1 >> log_tile_size) << 1;
cudaMalloc(&rank, sizeof(unsigned int) * index_len << 1);
unsigned int* split = rank + index_len;
unsigned int b = n + (tile_size << 1) - 1 >> log_tile_size + 1, s = tile_size << 1, m = n + tile_size - 1 >> log_tile_size;
binoticsort_gpu<<<b, s >> 1>>>(data, n);
for (b = b + 1 >> 1; s < n; s <<= 1, b = b + 1 >> 1)
{
unsigned int p = s + tile_size - 1 >> log_tile_size, k = (n - 1) / s + 1 >> 1;
rank_gpu<<<m + tile_size - 1 >> log_tile_size, tile_size>>>(data, n, s, rank);
ranksort_gpu<<<dim3(p + tile_size - 1 >> log_tile_size, std::min(k, 32768u), k + 32767 >> 15), tile_size>>>(rank, n, s, split);
if (!((n - 1) / s & 1))
++k;
mergesegment_gpu<<<dim3(p << 1, std::min(k, 32768u), k + 32767 >> 15), tile_size>>>(data, split, n, s, tmp);
std::swap(data, tmp);
}
cudaDeviceSynchronize();
cudaFree(tmp);
cudaFree(rank);
return data;
}投稿作者为『自动驾驶之心知识星球』特邀嘉宾，欢迎加入交流！重磅，自动驾驶之心科研论文辅导来啦，申博、CCF系列、SCI、EI、毕业论文、比赛辅导等多个方向，欢迎联系我们！
① 全网独家视频课程
BEV感知、BEV模型部署、BEV目标跟踪、毫米波雷达视觉融合、多传感器标定、多传感器融合、多模态3D目标检测、车道线检测、轨迹预测、在线高精地图、世界模型、点云3D目标检测、目标跟踪、Occupancy、cuda与TensorRT模型部署、大模型与自动驾驶、Nerf、语义分割、自动驾驶仿真、传感器部署、决策规划、轨迹预测等多个方向学习视频（扫码即可学习）
② 国内首个自动驾驶学习社区
国内最大最专业，近3000人的交流社区，已得到大多数自动驾驶公司的认可！涉及30+自动驾驶技术栈学习路线，从0到一带你入门自动驾驶感知（2D/3D检测、语义分割、车道线、BEV感知、Occupancy、多传感器融合、多传感器标定、目标跟踪）、自动驾驶定位建图（SLAM、高精地图、局部在线地图）、自动驾驶规划控制/轨迹预测等领域技术方案、大模型、端到端等，更有行业动态和岗位发布！欢迎扫描下方二维码，加入自动驾驶之心知识星球，这是一个真正有干货的地方，与领域大佬交流入门、学习、工作、跳槽上的各类难题，日常分享论文+代码+视频
③【自动驾驶之心】技术交流群
自动驾驶之心是首个自动驾驶开发者社区，聚焦感知、定位、融合、规控、标定、端到端、仿真、产品经理、自动驾驶开发、自动标注与数据闭环多个方向，目前近60+技术交流群，欢迎加入！扫码添加汽车人助理微信邀请入群，备注：学校/公司+方向+昵称（快速入群方式）
④【自动驾驶之心】全平台矩阵
扫一扫
741
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
