前言

前两章学了归约和归并，这章学排序——把无序数据变成有序数据。排序是最基础的算法之一，几乎所有领域都用得到。GPU 排序的挑战在于：传统排序算法（如快速排序）的递归和数据依赖不适合 GPU 的 SIMT 模型。第十三章介绍适合 GPU 的排序算法，重点是基数排序（Radix Sort）——利用前缀和实现高效并行排序。

📦 配套资源：本系列文章配有完整的 GitHub 仓库，包含每章的练习题解答、CUDA 代码实现和详细注释。所有代码都经过测试，可以直接运行。

排序算法概览

基于比较的排序

算法 平均复杂度 最坏复杂度 空间 稳定性
冒泡排序 O(n²) O(n²) O(1) 稳定
插入排序 O(n²) O(n²) O(1) 稳定
快速排序 O(n log n) O(n²) O(log n) 不稳定
归并排序 O(n log n) O(n log n) O(n) 稳定
堆排序 O(n log n) O(n log n) O(1) 不稳定

理论下界：基于比较的排序最优是 O(n log n)。

非比较排序

算法 复杂度 适用类型
计数排序 O(n + k) 小范围整数
基数排序 O(d(n + k)) 固定位数整数
桶排序 O(n + k) 均匀分布

突破下界：非比较排序可以达到 O(n)，但有类型限制。

GPU 排序的选择

快速排序不适合 GPU：

  • 递归深度不确定
  • 分区不平衡导致负载不均
  • 数据依赖导致分支发散

GPU 友好的算法：

  • 基数排序：规则的数据访问模式，利用前缀和
  • 归并排序：确定的步骤数，可并行归并
  • 双调排序：完全数据无关的比较网络

基数排序

核心思想

按位排序，从最低位到最高位：

1
2
3
4
5
6
7
8
9
10
11
12
原始:    [329, 457, 657, 839, 436, 720, 355]

按个位排序（第0位）:
         [720, 355, 436, 457, 657, 329, 839]

按十位排序（第1位）:
         [720, 329, 436, 839, 355, 457, 657]

按百位排序（第2位）:
         [329, 355, 436, 457, 657, 720, 839]

完成！

关键：每一轮排序必须是稳定的，保持相同键值元素的相对顺序。

位排序与计数排序

每一轮按某一位排序，使用计数排序：

1
2
3
1. 统计每个值的出现次数
2. 计算前缀和得到每个值的起始位置
3. 按原顺序分配到输出位置

对于二进制位（0 或 1）：

1
2
3
4
5
输入:    [1, 0, 1, 1, 0, 0, 1, 0]
计数:    count[0]=4, count[1]=4
前缀和:  pos[0]=0, pos[1]=4
分配:    输出位置 = pos[bit]++
结果:    [0, 0, 0, 0, 1, 1, 1, 1]

多位一起处理

一次处理多位（如 4 位）可以减少轮数：

1
2
3
4
32 位整数:
  1 位/轮: 32 轮
  4 位/轮: 8 轮
  8 位/轮: 4 轮

代价是每轮需要更大的计数数组（2^b 个桶）。

并行基数排序

单 Block 版本

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
#define RADIX_BITS 4
#define RADIX (1 << RADIX_BITS)  // 16
#define RADIX_MASK (RADIX - 1)

__global__ void radix_sort_block(unsigned int *data, int n, int bit) {
    __shared__ unsigned int s_data[BLOCK_SIZE];
    __shared__ int s_count[RADIX];
    __shared__ int s_offset[RADIX];
    
    int tid = threadIdx.x;
    
    // 加载数据
    s_data[tid] = (tid < n) ? data[tid] : 0xFFFFFFFF;
    __syncthreads();
    
    // 统计每个桶的计数（这里简化，实际需要原子或规约）
    if (tid < RADIX) s_count[tid] = 0;
    __syncthreads();
    
    int digit = (s_data[tid] >> bit) & RADIX_MASK;
    atomicAdd(&s_count[digit], 1);
    __syncthreads();
    
    // 前缀和得到偏移
    // ... exclusive scan on s_count → s_offset ...
    
    // 分配到输出位置
    int pos = atomicAdd(&s_offset[digit], 1);  // 动态分配
    
    __shared__ unsigned int s_temp[BLOCK_SIZE];
    s_temp[pos] = s_data[tid];
    __syncthreads();
    
    // 写回
    s_data[tid] = s_temp[tid];
    __syncthreads();
    
    if (tid < n) data[tid] = s_data[tid];
}

大规模基数排序

对于超过单 Block 的数据，需要分阶段：

阶段 1：局部直方图

每个 Block 统计自己区域的桶计数。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
__global__ void compute_histograms(unsigned int *data, int *histograms, 
                                    int n, int bit) {
    __shared__ int local_hist[RADIX];
    
    // 初始化
    if (threadIdx.x < RADIX) local_hist[threadIdx.x] = 0;
    __syncthreads();
    
    // 统计
    int i = blockIdx.x * blockDim.x + threadIdx.x;
    if (i < n) {
        int digit = (data[i] >> bit) & RADIX_MASK;
        atomicAdd(&local_hist[digit], 1);
    }
    __syncthreads();
    
    // 写入全局直方图
    if (threadIdx.x < RADIX) {
        histograms[blockIdx.x * RADIX + threadIdx.x] = local_hist[threadIdx.x];
    }
}

阶段 2：前缀和

对所有 Block 的直方图做全局前缀和，得到每个 Block 每个桶的全局起始位置。

1
Block 0 桶 k 的起始位置 = Σ(Block < 0, 桶 < k 的计数) + Σ(Block < 0, 桶 k 的计数)

这是一个二维前缀和问题。

阶段 3：重排（Scatter）

根据计算出的位置，把元素移动到正确位置。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
__global__ void scatter(unsigned int *in, unsigned int *out, 
                        int *prefix, int n, int bit) {
    __shared__ int local_prefix[RADIX];
    
    // 加载本 Block 的前缀
    if (threadIdx.x < RADIX) {
        local_prefix[threadIdx.x] = prefix[blockIdx.x * RADIX + threadIdx.x];
    }
    __syncthreads();
    
    int i = blockIdx.x * blockDim.x + threadIdx.x;
    if (i < n) {
        int digit = (in[i] >> bit) & RADIX_MASK;
        int pos = atomicAdd(&local_prefix[digit], 1);
        out[pos] = in[i];
    }
}

完整流程

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
void radix_sort(unsigned int *data, int n) {
    unsigned int *d_data, *d_temp;
    int *d_histograms, *d_prefix;
    
    int num_blocks = (n + BLOCK_SIZE - 1) / BLOCK_SIZE;
    
    cudaMalloc(&d_data, n * sizeof(unsigned int));
    cudaMalloc(&d_temp, n * sizeof(unsigned int));
    cudaMalloc(&d_histograms, num_blocks * RADIX * sizeof(int));
    cudaMalloc(&d_prefix, num_blocks * RADIX * sizeof(int));
    
    cudaMemcpy(d_data, data, n * sizeof(unsigned int), cudaMemcpyHostToDevice);
    
    for (int bit = 0; bit < 32; bit += RADIX_BITS) {
        // 阶段 1：计算直方图
        compute_histograms<<<num_blocks, BLOCK_SIZE>>>
            (d_data, d_histograms, n, bit);
        
        // 阶段 2：全局前缀和
        exclusive_scan(d_histograms, d_prefix, num_blocks * RADIX);
        
        // 阶段 3：重排
        scatter<<<num_blocks, BLOCK_SIZE>>>
            (d_data, d_temp, d_prefix, n, bit);
        
        // 交换缓冲区
        swap(d_data, d_temp);
    }
    
    cudaMemcpy(data, d_data, n * sizeof(unsigned int), cudaMemcpyDeviceToHost);
}

优化技术

避免原子操作

阶段 3 的 atomicAdd 是瓶颈。可以用本地排序 + 前缀和替代：

1
2
3
4
1. 每个线程处理多个元素
2. 本地计算每个桶的元素个数
3. Block 内前缀和得到每个线程的起始偏移
4. 无冲突地写入输出

向量化加载

每次加载 4 个元素：

1
2
uint4 data4 = *reinterpret_cast<uint4*>(&data[i]);
// 处理 data4.x, data4.y, data4.z, data4.w

减少内存事务数。

分桶优化

对于大范围数据，先按高位分成大桶，再分别排序小桶：

1
2
第一遍：按高 8 位分成 256 个桶
第二遍：对每个桶按剩余 24 位排序

减少每轮需要处理的数据量。

归并排序

GPU 归并排序流程

1
2
1. Block 内排序（用共享内存）
2. 跨 Block 归并（用并行归并）

Block 内排序可以用：

  • 比特序列排序（Bitonic Sort）
  • 奇偶归并排序
  • 直接归并排序

Block 内排序

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
__global__ void block_sort(int *data, int n) {
    __shared__ int s_data[BLOCK_SIZE];
    
    int tid = threadIdx.x;
    int gid = blockIdx.x * blockDim.x + tid;
    
    // 加载
    s_data[tid] = (gid < n) ? data[gid] : INT_MAX;
    __syncthreads();
    
    // 双调排序（Bitonic Sort）
    for (int k = 2; k <= BLOCK_SIZE; k *= 2) {
        for (int j = k / 2; j > 0; j /= 2) {
            int ixj = tid ^ j;
            if (ixj > tid) {
                bool ascending = ((tid & k) == 0);
                if ((s_data[tid] > s_data[ixj]) == ascending) {
                    // 交换
                    int temp = s_data[tid];
                    s_data[tid] = s_data[ixj];
                    s_data[ixj] = temp;
                }
            }
            __syncthreads();
        }
    }
    
    // 写回
    if (gid < n) data[gid] = s_data[tid];
}

跨 Block 归并

使用第十二章的并行归并技术：

1
2
3
4
5
6
7
8
9
10
11
12
13
void merge_sort_gpu(int *data, int n) {
    // 阶段 1：Block 内排序
    int num_blocks = (n + BLOCK_SIZE - 1) / BLOCK_SIZE;
    block_sort<<<num_blocks, BLOCK_SIZE>>>(data, n);
    
    // 阶段 2：逐层归并
    for (int width = BLOCK_SIZE; width < n; width *= 2) {
        int num_merges = (n + 2 * width - 1) / (2 * width);
        
        parallel_merge<<<num_merges, BLOCK_SIZE>>>(
            data, n, width);
    }
}

复杂度分析

阶段 工作量 并行度
Block 排序 O(n log B) n/B 个 Block
归并 k 层 O(n) × log(n/B) 递减
总计 O(n log n)

样本排序（Sample Sort）

思想

  1. 采样：从数据中随机选取样本
  2. 排序样本：得到分割点（splitters）
  3. 分区：根据分割点把数据分到不同桶
  4. 桶排序：每个桶独立排序
  5. 合并：按顺序拼接所有桶

优势

  • 负载均衡更好（采样保证桶大小接近）
  • 减少同步开销（桶间独立）

GPU 实现要点

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
void sample_sort(int *data, int n) {
    // 1. 采样
    int sample_size = min(n, 1024);
    sample_and_sort(data, n, samples, sample_size);
    
    // 2. 选择分割点
    int num_buckets = 256;
    select_splitters(samples, sample_size, splitters, num_buckets);
    
    // 3. 分区
    partition_to_buckets(data, n, splitters, bucket_ids, bucket_counts);
    
    // 4. 前缀和得到桶偏移
    exclusive_scan(bucket_counts, bucket_offsets, num_buckets);
    
    // 5. 重排到桶
    scatter_to_buckets(data, bucket_ids, bucket_offsets, temp, n);
    
    // 6. 桶内排序
    for (int b = 0; b < num_buckets; b++) {
        int start = bucket_offsets[b];
        int end = bucket_offsets[b + 1];
        sort_bucket<<<...>>>(temp + start, end - start);
    }
    
    // 7. 拷贝回
    cudaMemcpy(data, temp, n * sizeof(int), cudaMemcpyDeviceToDevice);
}

排序稳定性

为什么重要

稳定排序：相同键值的元素保持原有的相对顺序。

1
2
3
4
原始:    [(3,'a'), (1,'b'), (3,'c'), (2,'d')]
按数字稳定排序:
         [(1,'b'), (2,'d'), (3,'a'), (3,'c')]
                          ↑ 'a' 仍在 'c' 前面

应用：

  • 多关键字排序（先按次关键字，再按主关键字）
  • 数据库操作（保持插入顺序）
  • 基数排序的正确性依赖稳定性

GPU 排序的稳定性

算法 稳定性
基数排序 稳定
归并排序 稳定
双调排序 不稳定
快速排序 不稳定

基数排序天然稳定，是 GPU 排序的首选。

CUB/Thrust 库

使用 Thrust

1
2
3
4
5
6
7
8
9
10
11
12
13
#include <thrust/sort.h>
#include <thrust/device_vector.h>

void sortWithThrust(int *d_data, int n) {
    thrust::device_ptr<int> d_ptr(d_data);
    thrust::sort(d_ptr, d_ptr + n);
}

// 稳定排序
thrust::stable_sort(d_ptr, d_ptr + n);

// 按键排序
thrust::sort_by_key(d_keys, d_keys + n, d_values);

使用 CUB

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
#include <cub/cub.cuh>

void sortWithCub(unsigned int *d_data, int n) {
    void *d_temp = nullptr;
    size_t temp_bytes = 0;
    
    // 确定临时存储大小
    cub::DeviceRadixSort::SortKeys(d_temp, temp_bytes, 
                                    d_data, d_data, n);
    
    cudaMalloc(&d_temp, temp_bytes);
    
    // 执行排序
    cub::DeviceRadixSort::SortKeys(d_temp, temp_bytes,
                                    d_data, d_data, n);
    
    cudaFree(d_temp);
}

// 键值对排序
cub::DeviceRadixSort::SortPairs(d_temp, temp_bytes,
                                 d_keys, d_keys,
                                 d_values, d_values, n);

性能对比

以 1 亿个 32 位整数为例（RTX 3090）：

库/实现 时间 吞吐量
std::sort (CPU) ~8 秒 50 M/s
Thrust sort ~30 ms 3.3 G/s
CUB RadixSort ~15 ms 6.6 G/s

GPU 排序比 CPU 快 200-500 倍。

性能优化总结

优化层次

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
┌─────────────────────────────────────────────────┐
│ Level 4: 算法选择                                │
│   - 基数排序（整数最快）                         │
│   - 归并排序（通用稳定）                         │
│   - 样本排序（负载均衡）                         │
├─────────────────────────────────────────────────┤
│ Level 3: 数据局部性                              │
│   - Block 内共享内存排序                         │
│   - 向量化加载/存储                              │
├─────────────────────────────────────────────────┤
│ Level 2: 减少同步                                │
│   - Warp 内排序无需同步                          │
│   - 减少 Kernel 启动次数                         │
├─────────────────────────────────────────────────┤
│ Level 1: 利用前缀和                              │
│   - 高效的直方图和分配                           │
│   - 避免原子操作争用                             │
└─────────────────────────────────────────────────┘

选择指南

数据类型 推荐算法
32 位整数 基数排序
64 位整数 基数排序
浮点数 基数排序
结构体/自定义 归并排序
需要稳定性 基数/归并

浮点数需要特殊处理符号位和指数。

小结

第十三章系统讲解 GPU 排序：

基数排序：GPU 上最快的整数排序。按位处理，利用前缀和确定输出位置。多轮（32 位需要 8 轮 4 位）但每轮高度并行。

三阶段流程：直方图统计 → 前缀和计算 → 元素重排。每个阶段都是高效的并行操作。

归并排序：通用且稳定。Block 内用双调排序，跨 Block 用并行归并。适合任意可比较类型。

样本排序：通过采样实现负载均衡。适合数据分布不均匀的情况。

稳定性：基数排序天然稳定，多关键字排序和数据库操作必需。

库优先：CUB 的 RadixSort 高度优化，生产环境直接使用。理解原理有助于特殊需求的定制。

排序是并行计算的试金石，把前面学的归约、前缀和、归并等技术综合运用。掌握 GPU 排序，就掌握了这些基础原语的实战应用。

参考资料：

  • Hwu, W., Kirk, D., & El Hajj, I. (2022). Programming Massively Parallel Processors: A Hands-on Approach (4th Edition). Morgan Kaufmann.
  • Satish, N., Harris, M., & Garland, M. (2009). Designing efficient sorting algorithms for manycore GPUs. IPDPS.
  • NVIDIA CUB Library - Radix Sort

本文 GitHub 仓库: https://github.com/psmarter/PMPP-Learning