并行排序：从归并网络到 GPU 双调排序

单线程排序的优化已经被人类做到了极致。pdqsort 在随机数据上几乎跑满了单核的指令吞吐量，Timsort 在真实数据上利用预排序结构把比较次数压到了理论下界附近。但当你面对 10 亿条记录，单核即使每秒处理 2 亿次比较，也需要将近 2 分钟。而同一时间，你机器上的 63 个核心和那块有 10000 个 CUDA core 的 GPU 都在闲着。

并行排序的核心问题不是”怎么排得快”——而是”怎么把排序这个天然带有数据依赖的任务，拆成大量独立的并行单元”。这个问题的答案，从 1968 年 Batcher 的排序网络理论，到 2024 年 NVIDIA CUB 库的 DeviceRadixSort 实现，跨越了半个世纪的计算机科学。

本文从理论到工程，完整覆盖并行排序的各个层面。

一、比较器网络与 0-1 原理

1.1 什么是排序网络

排序网络（sorting network）是一种由比较器（comparator）组成的固定连接网络。每个比较器接收两根导线上的值，将较小的输出到上方导线，较大的输出到下方导线（或反之）。关键性质：网络的拓扑结构在编译时就确定了，不依赖于输入数据。

形式化定义：一个 n 输入的比较器网络是一系列比较器 (i, j) 的有序列表，其中 0 <= i < j < n。比较器 (i, j) 的语义是：

if a[i] > a[j]:
    swap(a[i], a[j])

如果对任意输入序列，网络输出都是排好序的，则称该网络为排序网络。

排序网络的两个核心度量： - 深度（depth）：将可并行执行的比较器分层后的层数，决定了并行时间复杂度。 - 大小（size）：比较器总数，决定了总工作量。

1.2 0-1 原理

0-1 原理（Zero-One Principle）是排序网络理论的基石：

定理：如果一个比较器网络能正确排序所有 0-1 序列（即每个元素只能是 0 或 1），则它能正确排序任意可比较元素的序列。

证明思路使用反证法。假设网络对某个非 0-1 输入 x 排序失败，即存在输出位置 i < j 使得 x’[i] > x’[j]。定义阈值函数 f(v) = (v >= x’[j]) ? 1 : 0。可以证明，对输入序列逐元素应用 f 得到的 0-1 序列，在同一网络中的行为与原序列完全一致——因为比较器对 f(a), f(b) 的结果与对 a, b 的结果通过 f 映射是一致的。但 f(x’[i]) = 1, f(x’[j]) = 0，所以这个 0-1 序列也排序失败，与假设矛盾。

0-1 原理的价值在于：验证一个 n 输入排序网络的正确性，只需要测试 2^n 个 0-1 输入，而不是所有 n! 种排列。虽然 2^n 仍然是指数级的，但在 n 较小时是可行的，而且更重要的是，它为排序网络的正确性证明提供了一个优雅的归纳框架。

1.3 排序网络的下界

AKS 网络（Ajtai-Komlós-Szemerédi, 1983）证明了存在深度为 O(log n)、大小为 O(n log n) 的排序网络，这在渐近意义上是最优的。但 AKS 网络的常数因子极大（实际中完全不可用）。

实际可用的排序网络的深度和大小：

对于实际大小（n <= 几千），Batcher 的两个网络是最佳选择。

二、Batcher 的奇偶归并网络

2.1 奇偶归并的递归结构

Kenneth Batcher 在 1968 年提出了两种经典排序网络。先看奇偶归并网络（Odd-Even Merge Network）。

核心思想：给定两个已排序的序列 A[0..n-1] 和 B[0..n-1]，如何用比较器网络将它们归并？

递归构造： 1. 将 A 和 B 各分成奇数位置和偶数位置的子序列。 2. 递归归并偶数位置的子序列和奇数位置的子序列。 3. 最后对相邻元素做一轮比较-交换。

def odd_even_merge(lo, hi, step):
    """生成奇偶归并网络的比较器列表。
    归并 a[lo], a[lo+step], a[lo+2*step], ..., a[hi] 这个子序列。
    """
    n = (hi - lo) // step + 1
    if n <= 1:
        return
    if n == 2:
        yield (lo, lo + step)
        return
    # 递归归并偶数位置
    yield from odd_even_merge(lo, hi, step * 2)
    # 递归归并奇数位置
    yield from odd_even_merge(lo + step, hi, step * 2)
    # 最终的比较-交换
    i = lo + step
    while i + step <= hi:
        yield (i, i + step)
        i += step * 2


def odd_even_merge_sort(lo, hi):
    """生成奇偶归并排序网络。"""
    if hi - lo < 1:
        return
    mid = lo + (hi - lo) // 2
    yield from odd_even_merge_sort(lo, mid)
    yield from odd_even_merge_sort(mid + 1, hi)
    yield from odd_even_merge(lo, hi, 1)

2.2 复杂度分析

设 M(n) 为归并两个长度为 n/2 的已排序序列的比较器数。递归关系：

M(n) = 2 * M(n/2) + n/2 - 1
M(2) = 1

解为 M(n) = (n/2) * log_2(n/2) + 1 - n/2 ≈ (n/2) * log(n)。

整个排序的比较器数 S(n)：

S(n) = 2 * S(n/2) + M(n)

解为 S(n) = O(n * log^2(n))。深度为 O(log^2 n)。

2.3 与普通归并排序的区别

普通归并排序的归并操作是顺序的——你必须从两个序列的头部开始，逐个比较取小的那个。这个过程有严格的数据依赖：第 k 个输出的选择依赖于前 k-1 个选择。

奇偶归并网络的归并是”无条件的”——所有比较器的位置在编译时就确定了，不依赖数据。同一层的比较器可以完全并行执行。这正是排序网络能映射到硬件和 GPU 的关键。

三、双调排序（Bitonic Sort）

3.1 双调序列的定义

双调序列（bitonic sequence）是先单调递增再单调递减（或反之）的序列，或者可以通过循环移位变成先增后减的序列。形式化地说，序列 a[0..n-1] 是双调的，如果存在某个 k 使得：

a[0] <= a[1] <= ... <= a[k] >= a[k+1] >= ... >= a[n-1]

或者上述条件对 a 的某个循环移位成立。

例如：[1, 3, 5, 7, 6, 4, 2, 0] 是双调的（先增后减）。[5, 7, 6, 4, 2, 0, 1, 3] 也是双调的（循环移位后先增后减）。

3.2 双调分裂（Bitonic Split）

双调排序的核心操作是双调分裂：给定一个长度为 n 的双调序列，将每对 (a[i], a[i + n/2]) 做比较-交换，可以得到两个性质：

1. 前半部分的所有元素 <= 后半部分的所有元素。
2. 前半部分和后半部分各自仍然是双调序列。

这意味着我们可以递归地对两半分别做双调分裂，直到子序列长度为 1——此时整个序列就排好了。

双调分裂 n=8:
    对 i in 0..3:
        compare_swap(a[i], a[i+4])
    -> 两个长度为 4 的双调序列

双调分裂 n=4 (对两半分别):
    对 i in 0..1:
        compare_swap(a[i], a[i+2])
    -> 四个长度为 2 的双调序列

双调分裂 n=2 (对四段分别):
    compare_swap(a[i], a[i+1])
    -> 排好了

3.3 从双调归并到双调排序

双调归并（Bitonic Merge）：给定一个双调序列，通过 log_2(n) 步双调分裂将其排序。

双调排序的递归构造： 1. 将输入分成两半。 2. 递归地对左半部分按升序排序，右半部分按降序排序。 3. 此时整个序列是双调的（左半升序 + 右半降序）。 4. 对整个双调序列做双调归并。

/* 双调排序的比较-交换单元 */
void compare_swap(int *a, int i, int j, int dir) {
    /* dir=1 升序, dir=0 降序 */
    if (dir == (a[i] > a[j])) {
        int tmp = a[i];
        a[i] = a[j];
        a[j] = tmp;
    }
}

/* 双调归并：将双调序列排序 */
void bitonic_merge(int *a, int lo, int cnt, int dir) {
    if (cnt <= 1) return;
    int k = cnt / 2;
    for (int i = lo; i < lo + k; i++)
        compare_swap(a, i, i + k, dir);
    bitonic_merge(a, lo, k, dir);
    bitonic_merge(a, lo + k, k, dir);
}

/* 双调排序：递归构建双调序列，然后归并 */
void bitonic_sort(int *a, int lo, int cnt, int dir) {
    if (cnt <= 1) return;
    int k = cnt / 2;
    bitonic_sort(a, lo, k, 1);       /* 左半升序 */
    bitonic_sort(a, lo + k, k, 0);   /* 右半降序 */
    bitonic_merge(a, lo, cnt, dir);   /* 归并 */
}

3.4 复杂度

对于 n = 2^p 个元素： - 深度：p * (p + 1) / 2 = O(log^2 n) - 比较器数：n * p * (p + 1) / 4 = O(n * log^2 n) - 对 n = 2^20 ≈ 100 万，深度 = 210，比较器数 ≈ 1.1 亿

比较一下：普通归并排序的比较次数约 n * log_2(n) = 2000 万。双调排序做了大约 5 倍多的比较。但如果我们有足够的并行度（n/2 个处理器），双调排序只需要 210 步，而串行归并排序需要 2000 万步。这就是并行排序的本质交换：用更多的总工作量换取更少的关键路径长度。

3.5 非 2 的幂次大小的处理

实际应用中，数组大小不一定是 2 的幂。处理方法有几种：

1. 填充法：用极大值（正无穷）将数组填充到 2 的幂次大小。排序后忽略填充值。简单但浪费内存和计算。

2. 截断法：修改比较器逻辑，当索引越界时跳过比较。需要在 kernel 中加边界检查，但不浪费内存。

3. 分块法：将数组分成 2 的幂次大小的块分别排序，最后归并。工程上最常用。

/* 安全比较-交换，处理越界 */
void safe_compare_swap(int *a, int n, int i, int j, int dir) {
    if (j >= n) return;  /* 索引越界，跳过 */
    if (dir == (a[i] > a[j])) {
        int tmp = a[i];
        a[i] = a[j];
        a[j] = tmp;
    }
}

四、双调排序的 GPU 映射

4.1 为什么双调排序特别适合 GPU

GPU 的执行模型是 SIMT（Single Instruction, Multiple Threads）——同一个 warp 内的 32 个线程执行相同的指令。这对排序算法有几个关键约束：

1. 分支散度（warp divergence）惩罚巨大。如果同一个 warp 中的线程走不同的分支，两条路径都要执行——性能直接减半甚至更差。

2. 全局内存访问要合并（coalesced）。连续的线程访问连续的内存地址，可以合并成一次宽的内存事务；随机访问则产生大量独立事务。

3. 没有高效的跨线程通信。不同 block 之间只能通过全局内存同步。

双调排序完美匹配这些约束：

• 没有数据依赖的分支。每一步中，每个线程根据自己的线程 ID 就能计算出应该比较哪两个位置，比较方向也是确定的。同一个 warp 内的线程执行完全相同的指令序列——没有分支散度。

• 内存访问模式规则。在双调分裂的每一步中，线程 i 访问 a[i] 和 a[i XOR step]（其中 step 是当前步长）。这产生了一种蝴蝶式（butterfly pattern）的内存访问，步长固定且可预测。

• 同步点少且规则。每一步之后需要一次全局栅栏同步，然后进入下一步。没有复杂的任务调度或工作窃取。

4.2 完整的 CUDA 双调排序实现

#include <cuda_runtime.h>
#include <stdio.h>

/* ---- 全局内存版本的双调排序 kernel ---- */
__global__ void bitonic_sort_step(int *dev_values, int j, int k) {
    unsigned int i = threadIdx.x + blockDim.x * blockIdx.x;
    unsigned int ixj = i ^ j;  /* 蝴蝶式配对 */

    /* 每对只处理一次：i < ixj */
    if (ixj > i) {
        /* 判断排序方向：(i & k) == 0 -> 升序 */
        if ((i & k) == 0) {
            if (dev_values[i] > dev_values[ixj]) {
                int tmp = dev_values[i];
                dev_values[i] = dev_values[ixj];
                dev_values[ixj] = tmp;
            }
        } else {
            if (dev_values[i] < dev_values[ixj]) {
                int tmp = dev_values[i];
                dev_values[i] = dev_values[ixj];
                dev_values[ixj] = tmp;
            }
        }
    }
}

/* ---- 共享内存优化版本 ---- */
#define SHARED_SIZE_LIMIT 1024

__global__ void bitonic_sort_shared(int *dev_values, int n) {
    __shared__ int shared[SHARED_SIZE_LIMIT];

    unsigned int tid = threadIdx.x;
    unsigned int gid = blockIdx.x * blockDim.x + tid;

    /* 加载到共享内存 */
    shared[tid] = (gid < n) ? dev_values[gid] : INT_MAX;
    __syncthreads();

    /* 双调排序的所有阶段在共享内存中完成 */
    for (unsigned int k = 2; k <= blockDim.x; k <<= 1) {
        for (unsigned int j = k >> 1; j > 0; j >>= 1) {
            unsigned int ixj = tid ^ j;
            if (ixj > tid) {
                if ((tid & k) == 0) {
                    if (shared[tid] > shared[ixj]) {
                        int tmp = shared[tid];
                        shared[tid] = shared[ixj];
                        shared[ixj] = tmp;
                    }
                } else {
                    if (shared[tid] < shared[ixj]) {
                        int tmp = shared[tid];
                        shared[tid] = shared[ixj];
                        shared[ixj] = tmp;
                    }
                }
            }
            __syncthreads();
        }
    }

    /* 写回全局内存 */
    if (gid < n) {
        dev_values[gid] = shared[tid];
    }
}

/* ---- 混合策略：小块用共享内存，大步长用全局内存 ---- */
void bitonic_sort_gpu(int *h_values, int n) {
    int *d_values;
    size_t size = n * sizeof(int);

    cudaMalloc(&d_values, size);
    cudaMemcpy(d_values, h_values, size, cudaMemcpyHostToDevice);

    dim3 blocks(n / SHARED_SIZE_LIMIT);
    dim3 threads(SHARED_SIZE_LIMIT);

    /* 第一阶段：块内排序，使用共享内存 */
    bitonic_sort_shared<<<blocks, threads>>>(d_values, n);

    /* 第二阶段：跨块归并，使用全局内存 */
    for (int k = SHARED_SIZE_LIMIT * 2; k <= n; k <<= 1) {
        for (int j = k >> 1; j > 0; j >>= 1) {
            if (j >= SHARED_SIZE_LIMIT) {
                /* 大步长：全局内存 */
                bitonic_sort_step<<<blocks, threads>>>(d_values, j, k);
            } else {
                /* 小步长：可以用共享内存优化
                   这里简化处理，仍用全局内存 */
                bitonic_sort_step<<<blocks, threads>>>(d_values, j, k);
            }
        }
    }

    cudaMemcpy(h_values, d_values, size, cudaMemcpyDeviceToHost);
    cudaFree(d_values);
}

/* ---- 测试入口 ---- */
int main() {
    const int N = 1 << 20;  /* 100 万元素 */
    int *data = (int *)malloc(N * sizeof(int));

    srand(42);
    for (int i = 0; i < N; i++)
        data[i] = rand();

    bitonic_sort_gpu(data, N);

    /* 验证 */
    for (int i = 1; i < N; i++) {
        if (data[i - 1] > data[i]) {
            printf("排序错误 at index %d\n", i);
            free(data);
            return 1;
        }
    }
    printf("排序正确，%d 个元素\n", N);

    free(data);
    return 0;
}

4.3 性能分析与优化要点

双调排序在 GPU 上的性能瓶颈：

1. 全局内存带宽。每一步的双调分裂都需要读写所有 n 个元素。对 log^2(n) 步来说，总内存访问量是 O(n * log^2(n))。对 n = 100 万，这大约是 210 次全局内存扫描。GPU 的全局内存带宽虽然高（例如 A100 的 2TB/s），但内存访问仍然是主要瓶颈。

2. 共享内存优化。当步长 j 足够小（j < block 大小）时，一个 block 内的所有比较-交换可以在共享内存中完成，避免全局内存访问。这将一个 block 内的 log(block_size) 步合并成一次全局内存读 + 多次共享内存操作 + 一次全局内存写。

3. 寄存器优化。更激进地，当步长极小时（j < warp 大小），可以利用 warp shuffle 指令在寄存器级别交换数据，完全不经过共享内存。

/* 使用 warp shuffle 进行 warp 级别的双调排序 */
__device__ int warp_bitonic_sort(int val, int lane_id) {
    for (int k = 2; k <= 32; k <<= 1) {
        for (int j = k >> 1; j > 0; j >>= 1) {
            int partner = lane_id ^ j;
            int partner_val = __shfl_xor_sync(0xFFFFFFFF, val, j);
            bool ascending = ((lane_id & k) == 0);
            if (ascending) {
                val = (val > partner_val && partner > lane_id) ?
                      partner_val : val;
                val = (val < partner_val && partner < lane_id) ?
                      partner_val : val;
            } else {
                val = (val < partner_val && partner > lane_id) ?
                      partner_val : val;
                val = (val > partner_val && partner < lane_id) ?
                      partner_val : val;
            }
        }
    }
    return val;
}

4.4 双调排序 vs 基数排序在 GPU 上的对比

结论很明显：对于整数和浮点数排序，GPU 基数排序碾压双调排序。但对于自定义比较器或键值对排序，双调排序仍有其生态位。

五、并行归并排序与 Merge Path 算法

5.1 并行归并的挑战

传统归并排序的递归结构天然适合并行化——每一层的多个归并可以独立进行。但问题出在最后一层：当你需要归并两个各 n/2 大小的已排序序列时，传统的双指针归并是严格顺序的。

如何并行化单次归并操作？这是并行归并排序的核心问题。

5.2 Merge Path 算法

Merge Path（由 Odeh、Green、Mwassi、Shmueli 和 Birk 于 2012 年提出）是一种优雅的方法，可以将两个已排序序列的归并操作均匀地分配给 p 个处理器。

核心观察：考虑两个已排序序列 A[0..m-1] 和 B[0..n-1]。定义一个 (m+1) x (n+1) 的矩阵，其中 M[i][j] = (A[i-1] > B[j-1])。归并过程可以看作在这个矩阵中从 (0,0) 走到 (m,n) 的路径——每次向右走表示从 B 取一个元素，向下走表示从 A 取一个元素。

关键性质：这条路径是单调的（只能向右或向下走），且路径上方全是 1，下方全是 0。因此，对于任意对角线（i + j = k），路径与该对角线恰好交于一点，可以用二分查找在 O(log(min(m,n))) 时间内找到。

/* Merge Path 分区：找到输出位置 diag 对应的 A 和 B 的切分点 */
int merge_path(const int *A, int m, const int *B, int n, int diag) {
    int lo = (diag > n) ? diag - n : 0;
    int hi = (diag < m) ? diag : m;

    while (lo < hi) {
        int mid = lo + (hi - lo) / 2;
        /* A[mid] 与 B[diag - mid - 1] 比较 */
        if (A[mid] <= B[diag - mid - 1])
            lo = mid + 1;
        else
            hi = mid;
    }
    return lo;  /* A 中取 lo 个，B 中取 diag - lo 个 */
}

/* 并行归并：p 个线程各自处理 (m+n)/p 个输出元素 */
void parallel_merge(const int *A, int m, const int *B, int n,
                    int *out, int p) {
    int total = m + n;
    int chunk = (total + p - 1) / p;

    #pragma omp parallel for
    for (int tid = 0; tid < p; tid++) {
        int diag_start = tid * chunk;
        int diag_end = (tid + 1) * chunk;
        if (diag_end > total) diag_end = total;

        /* 找到这个线程负责的 A 和 B 的区间 */
        int a_start = merge_path(A, m, B, n, diag_start);
        int b_start = diag_start - a_start;
        int a_end = merge_path(A, m, B, n, diag_end);
        int b_end = diag_end - a_end;

        /* 顺序归并这个小区间 */
        int ai = a_start, bi = b_start, oi = diag_start;
        while (ai < a_end && bi < b_end) {
            if (A[ai] <= B[bi])
                out[oi++] = A[ai++];
            else
                out[oi++] = B[bi++];
        }
        while (ai < a_end) out[oi++] = A[ai++];
        while (bi < b_end) out[oi++] = B[bi++];
    }
}

5.3 work-efficient 分析

Merge Path 的工作量分析： - 每个线程的分区查找：O(log(min(m,n))) - 每个线程的归并：O((m+n)/p) - 总工作量：O(m+n + p * log(min(m,n))) - 并行时间：O((m+n)/p + log(min(m,n)))

当 p << (m+n) / log(min(m,n)) 时，这是 work-efficient 的——总工作量与串行归并的 O(m+n) 相差一个 log 因子。

5.4 GPU 上的 Merge Path

ModernGPU 库将 Merge Path 应用到了 GPU 上。关键优化：

1. 两级分区。先用 block 级别的 Merge Path 将大归并分成每个 block 的子任务（每个 block 处理几千个元素），然后在 block 内用共享内存完成归并。

2. 向量化加载。使用 128 位加载指令一次读 4 个 int，减少内存事务数。

3. 串行归并在共享内存中完成。每个线程在共享内存中归并一小段，避免全局内存的延迟。

六、多核 CPU 上的并行排序

6.1 GNU parallel mode

GCC 的 libstdc++ 提供了 parallel mode，通过 OpenMP 实现并行排序：

#include <parallel/algorithm>

void sort_parallel(std::vector<int>& v) {
    __gnu_parallel::sort(v.begin(), v.end());
}

/* 编译：g++ -fopenmp -D_GLIBCXX_PARALLEL sort.cpp */

GNU parallel mode 的排序策略：

• 数据量较小：直接调用串行排序（通常是 introsort）。
• 数据量中等（几十万到几百万）：并行多路归并排序。将数据分成 p 份（p 是核心数），每个核心独立排序，最后用多路归并合并。
• 数据量较大（千万以上）：并行采样排序。通过采样确定分区点，将数据分成 p 个桶，每个核心排序一个桶。

6.2 Intel TBB parallel_sort

Intel Threading Building Blocks (TBB) 的 tbb::parallel_sort 实现了并行快速排序：

#include <tbb/parallel_sort.h>
#include <vector>

void sort_with_tbb(std::vector<int>& v) {
    tbb::parallel_sort(v.begin(), v.end());
}

TBB parallel_sort 的实现策略：

1. 在顶层使用并行快速排序——每次 partition 后，将两个子问题作为 TBB task 提交。
2. 当子问题大小低于阈值（通常 500-2000 个元素）时，回退到串行排序。
3. 使用工作窃取（work stealing）调度器平衡负载——如果某个核心的 partition 不均匀导致一侧特别大，其他空闲核心可以窃取任务。

/* TBB 并行快排的概念性实现 */
template<typename Iter, typename Cmp>
void parallel_quicksort(Iter begin, Iter end, Cmp cmp, int depth) {
    if (end - begin < 2000 || depth > 16) {
        std::sort(begin, end, cmp);
        return;
    }

    auto pivot = median_of_three(begin, end, cmp);
    auto mid = partition(begin, end, pivot, cmp);

    tbb::parallel_invoke(
        [&]{ parallel_quicksort(begin, mid, cmp, depth + 1); },
        [&]{ parallel_quicksort(mid, end, cmp, depth + 1); }
    );
}

6.3 OpenMP 并行归并排序的完整实现

#include <stdlib.h>
#include <string.h>
#include <omp.h>

/* 串行归并 */
static void merge(int *arr, int *tmp, int lo, int mid, int hi) {
    int i = lo, j = mid, k = lo;
    while (i < mid && j < hi) {
        if (arr[i] <= arr[j])
            tmp[k++] = arr[i++];
        else
            tmp[k++] = arr[j++];
    }
    while (i < mid) tmp[k++] = arr[i++];
    while (j < hi)  tmp[k++] = arr[j++];
    memcpy(arr + lo, tmp + lo, (hi - lo) * sizeof(int));
}

/* 并行归并排序 */
void parallel_merge_sort(int *arr, int *tmp, int lo, int hi, int depth) {
    if (hi - lo < 2) return;
    if (hi - lo < 4096 || depth > 4) {
        /* 小段或递归过深：回退串行排序 */
        /* 使用插入排序或标准库排序 */
        for (int i = lo + 1; i < hi; i++) {
            int key = arr[i];
            int j = i - 1;
            while (j >= lo && arr[j] > key) {
                arr[j + 1] = arr[j];
                j--;
            }
            arr[j + 1] = key;
        }
        return;
    }

    int mid = lo + (hi - lo) / 2;

    #pragma omp task shared(arr, tmp) if(depth < 4)
    parallel_merge_sort(arr, tmp, lo, mid, depth + 1);

    #pragma omp task shared(arr, tmp) if(depth < 4)
    parallel_merge_sort(arr, tmp, mid, hi, depth + 1);

    #pragma omp taskwait
    merge(arr, tmp, lo, mid, hi);
}

void sort(int *arr, int n) {
    int *tmp = (int *)malloc(n * sizeof(int));
    #pragma omp parallel
    {
        #pragma omp single
        parallel_merge_sort(arr, tmp, 0, n, 0);
    }
    free(tmp);
}

6.4 负载均衡问题

并行快速排序的最大问题是负载不均衡。快速排序的 partition 将数组分成两部分，比例取决于 pivot 的选择。即使用三数取中，worst-case 的不均衡也是 1:n-2。这意味着：

• 如果第一次 partition 极度不均衡（比如 1:99），一个核心分到 99% 的数据，其他核心很快完成后只能空等。
• 工作窃取可以缓解但不能根治——窃取有开销，而且如果不均衡发生在递归的顶层，影响最大。

并行归并排序没有这个问题——每次递归都是精确的二分。但归并排序需要额外的 O(n) 空间，且缓存行为不如快速排序。

七、采样排序（Sample Sort）

7.1 分布式环境的需求

在多机分布式环境中（MPI 集群、MapReduce），数据分布在不同机器上，通信成本远高于计算成本。并行快速排序和并行归并排序都需要大量的数据移动——partition 时要在机器间交换元素，归并时要从多台机器拉取数据。

采样排序的核心思想：通过采样估算数据分布，选择合适的分区点（splitters），使得每个处理器分配到大致相等数量的元素。一次通信确定分区方案，一次全对全通信（all-to-all）完成数据重分布，然后每个处理器独立排序本地数据。

7.2 算法流程

给定 p 个处理器，每个持有 n/p 个元素：

1. 每个处理器对本地数据排序
2. 每个处理器从本地选取 s 个等间距样本（oversampling factor s，通常 s = log(n) 或 s = p）
3. 收集所有 p*s 个样本到一个处理器
4. 对样本排序，选取 p-1 个分割点（splitters）
5. 广播 splitters 到所有处理器
6. 每个处理器根据 splitters 将本地数据分成 p 个桶
7. All-to-all 通信：处理器 i 的第 j 个桶发送给处理器 j
8. 每个处理器对收到的数据排序

#include <mpi.h>
#include <algorithm>
#include <vector>

void sample_sort(std::vector<int>& local_data, MPI_Comm comm) {
    int rank, nprocs;
    MPI_Comm_rank(comm, &rank);
    MPI_Comm_size(comm, &nprocs);

    /* 步骤 1：本地排序 */
    std::sort(local_data.begin(), local_data.end());
    int local_n = local_data.size();

    /* 步骤 2：选取本地样本 */
    int s = nprocs;  /* oversampling factor */
    std::vector<int> local_samples(s);
    for (int i = 0; i < s; i++)
        local_samples[i] = local_data[i * local_n / s];

    /* 步骤 3：收集所有样本到 rank 0 */
    std::vector<int> all_samples;
    if (rank == 0) all_samples.resize(nprocs * s);
    MPI_Gather(local_samples.data(), s, MPI_INT,
               all_samples.data(), s, MPI_INT, 0, comm);

    /* 步骤 4：选取 splitters */
    std::vector<int> splitters(nprocs - 1);
    if (rank == 0) {
        std::sort(all_samples.begin(), all_samples.end());
        for (int i = 0; i < nprocs - 1; i++)
            splitters[i] = all_samples[(i + 1) * s];
    }

    /* 步骤 5：广播 splitters */
    MPI_Bcast(splitters.data(), nprocs - 1, MPI_INT, 0, comm);

    /* 步骤 6：根据 splitters 分桶 */
    std::vector<std::vector<int>> buckets(nprocs);
    for (int v : local_data) {
        int dest = std::lower_bound(splitters.begin(),
                                     splitters.end(), v)
                   - splitters.begin();
        buckets[dest].push_back(v);
    }

    /* 步骤 7：All-to-all 通信 */
    std::vector<int> send_counts(nprocs), recv_counts(nprocs);
    for (int i = 0; i < nprocs; i++)
        send_counts[i] = buckets[i].size();
    MPI_Alltoall(send_counts.data(), 1, MPI_INT,
                 recv_counts.data(), 1, MPI_INT, comm);

    std::vector<int> send_displs(nprocs), recv_displs(nprocs);
    send_displs[0] = recv_displs[0] = 0;
    for (int i = 1; i < nprocs; i++) {
        send_displs[i] = send_displs[i-1] + send_counts[i-1];
        recv_displs[i] = recv_displs[i-1] + recv_counts[i-1];
    }

    int total_recv = recv_displs[nprocs-1] + recv_counts[nprocs-1];
    std::vector<int> send_buf, recv_buf(total_recv);
    for (auto& b : buckets)
        send_buf.insert(send_buf.end(), b.begin(), b.end());

    MPI_Alltoallv(send_buf.data(), send_counts.data(),
                  send_displs.data(), MPI_INT,
                  recv_buf.data(), recv_counts.data(),
                  recv_displs.data(), MPI_INT, comm);

    /* 步骤 8：本地排序 */
    std::sort(recv_buf.begin(), recv_buf.end());
    local_data = std::move(recv_buf);
}

7.3 关于过采样因子的选择

过采样因子 s 的选择是一个精妙的权衡：

• s 太小：splitters 不够精确，导致桶大小不均衡。极端情况下某个桶可能分到 80% 的数据。
• s 太大：采样和排序样本的开销增加，通信量也增大。

理论结果：当 s = O(log(n) / epsilon^2) 时，可以保证每个桶的大小在 n/p * (1 +/- epsilon) 范围内，概率大于 1 - 1/n。

实践中，s = max(p, 64) 通常就足够了。

7.4 与 MapReduce 的关系

Hadoop 的 TeraSort（赢得 2008 年排序基准测试的实现）本质上就是采样排序：

1. Map 阶段：采样估算数据分布，确定分区点。
2. Shuffle 阶段：根据分区点将数据路由到对应的 Reducer（即 all-to-all 通信）。
3. Reduce 阶段：每个 Reducer 对本地数据排序。

Spark 的 sortByKey 也使用类似的采样分区策略。

八、GPU 排序的工业实现

8.1 CUB DeviceRadixSort

NVIDIA CUB 库的 DeviceRadixSort 是目前 GPU 上最快的整数排序实现之一。它的核心思想：

1. 按位分解。对 32 位整数，分成 8 趟（每趟处理 4 位，即 16 个桶）或 4 趟（每趟 8 位，256 个桶）。

2. 每趟是一次前缀和 + 散射。

   • 统计每个桶中元素的数量（直方图）。
   • 对桶大小做前缀和，得到每个桶的起始位置。
   • 将元素散射到正确位置。

3. 优化细节：

   • Warp 级别的直方图，减少原子操作冲突。
   • 每个 block 先做本地的基数排序（使用共享内存），然后才写回全局内存。
   • 对最后几位可以用计数排序一步完成。

#include <cub/cub.cuh>

void sort_with_cub(int *d_keys, int n) {
    void *d_temp = nullptr;
    size_t temp_bytes = 0;

    /* 第一次调用确定临时空间大小 */
    cub::DeviceRadixSort::SortKeys(d_temp, temp_bytes, d_keys,
                                    d_keys, n);
    cudaMalloc(&d_temp, temp_bytes);

    /* 第二次调用执行排序 */
    cub::DeviceRadixSort::SortKeys(d_temp, temp_bytes, d_keys,
                                    d_keys, n);
    cudaFree(d_temp);
}

8.2 Thrust sort

Thrust 是 NVIDIA 的高级 C++ 并行算法库，提供了类 STL 的接口：

#include <thrust/sort.h>
#include <thrust/device_vector.h>

void sort_with_thrust(thrust::device_vector<int>& d_vec) {
    thrust::sort(d_vec.begin(), d_vec.end());
}

Thrust 的 sort 对整数类型内部调用 CUB 的基数排序，对自定义类型使用归并排序。

8.3 ModernGPU 的归并排序

ModernGPU（由 Sean Baxter 编写）实现了高效的 GPU 归并排序，其核心是前面讨论的 Merge Path 算法。

实现层次： 1. Warp 级别：每个线程持有若干元素，warp 内用寄存器级别的排序网络排序。 2. Block 级别：block 内的线程通过共享内存归并各自的有序段。 3. Device 级别：block 间通过全局内存归并，使用 Merge Path 做均匀分区。

ModernGPU 的优势在于它对任意比较器都能高效工作，不像基数排序那样局限于整数类型。

8.4 GPU 排序的选择指南

输入是整数/浮点数？
  ├─ 是 → CUB DeviceRadixSort（最快）
  └─ 否 → 元素大小？
           ├─ 小（key-value pair < 16 字节）→ Thrust sort
           └─ 大（复杂结构体）→ 先排序索引，再 gather

九、实测对比

9.1 测试环境

以下数据基于典型硬件配置，用于说明数量级关系：

• CPU：AMD Ryzen 9 7950X（16 核 32 线程，5.7GHz boost）
• GPU：NVIDIA RTX 4090（16384 CUDA cores，1TB/s 内存带宽）
• 内存：64GB DDR5-5600
• 数据：随机 32 位整数

9.2 性能数据

注意：GPU 时间不包含 CPU-GPU 数据传输。如果包含传输时间，n=10^6 时 GPU 排序通常比 CPU 更慢（传输延迟约 0.5ms + 带宽瓶颈）。

9.3 加速比分析

n = 10^8 时的加速比（相对于 std::sort）：

pdqsort:              1.4x
GNU parallel (16核):  13.5x
TBB parallel (16核):  15.6x
Thrust (GPU):         300x
CUB radix (GPU):      525x

几个观察：

1. 多核 CPU 的加速比接近但达不到核心数（16）。原因是归并阶段的内存带宽瓶颈——16 个核心共享内存控制器，排序的最后几次归并变成了内存带宽受限。

2. GPU 的加速比远超核心数比。这是因为 GPU 不仅核心多，而且内存带宽高出一个数量级（1TB/s vs 80GB/s）。对于排序这种内存密集型任务，带宽优势被充分利用。

3. 当 n 足够大时（10^9），GPU 的优势非常明显。当 n 较小时（10^6），CPU-GPU 传输开销可能抵消计算优势。

9.4 含传输时间的真实场景

在实际应用中，数据通常在 CPU 内存中。需要考虑 PCIe 传输：

PCIe 4.0 x16 带宽：约 25GB/s
10^8 个 int = 400MB
传输时间（双向）：约 32ms

含传输的 GPU 排序时间：32ms + 16ms = 48ms
CPU 16核并行排序：620ms

即使含传输，GPU 仍快 12 倍以上。
但对 n = 10^6（4MB）：
传输时间：约 0.3ms
GPU 排序：0.3ms + 0.8ms = 1.1ms
CPU 16核并行：8ms

仍然是 GPU 更快，但差距缩小到 7 倍。

十、并行排序的 Scalability 分析

10.1 Amdahl 定律

Amdahl 定律指出，如果程序中有比例为 f 的部分必须串行执行，则使用 p 个处理器的最大加速比为：

S(p) = 1 / (f + (1 - f) / p)

当 p -> infinity 时，S -> 1/f。即串行部分决定了加速比的上限。

10.2 排序中的串行瓶颈

并行排序中的串行部分包括：

1. 递归树的顶层。并行归并排序/快速排序的递归树，前 log_2(p) 层才能利用所有核心。第一次 partition 或第一次归并只用 1-2 个核心。

2. 采样和分区决策。采样排序中，收集样本、排序样本、选择分割点的过程通常在单个处理器上完成。

3. 结果合并。最后一次归并，或最后的数据收集。

对于 n 个元素、p 个处理器的并行归并排序，串行部分大约是 O(n) 的最终归并（如果不用 Merge Path 并行化的话），占总工作量 O(n log n) 的比例 f = 1/log(n)。

f = 1 / log_2(n)
对 n = 10^9, f ≈ 1/30 = 3.3%
S_max = 1 / 0.033 = 30

即使有无限核心，这种实现的加速比也不超过 30 倍。

10.3 Gustafson 定律与弱扩展

Amdahl 定律假设问题规模固定。Gustafson 定律考虑另一种场景：增加处理器时，也增加问题规模（弱扩展，weak scaling）。

S_gustafson(p) = p - f * (p - 1) ≈ p * (1 - f)

对于排序，如果 n 随 p 线性增长（每个处理器始终处理相同数量的数据），加速比可以近似线性增长。这是分布式排序（如采样排序）的真实场景——你增加机器不是为了更快地排同样多的数据，而是为了排更多的数据。

10.4 内存带宽瓶颈

在多核 CPU 上，排序的 scalability 通常不是被计算瓶颈限制，而是被内存带宽限制。

分析：排序 n 个 4 字节整数，总数据量 4n 字节。归并排序需要约 4n * log_2(n) 字节的内存访问（每一层读写所有数据）。

n = 10^8, 4 字节整数
内存访问量 ≈ 400MB * 27 ≈ 10.8GB
DDR5-5600 带宽 ≈ 80GB/s
理论最短时间 ≈ 135ms

实测 16 核并行：~600ms
效率 ≈ 22.5%

效率不高的原因：不是所有内存访问都能打满带宽——cache miss 率、TLB miss、NUMA 跨节点访问都会降低有效带宽。

10.5 NUMA 的影响

在多 socket 系统上，NUMA（Non-Uniform Memory Access）效应显著。跨 NUMA 节点的内存访问延迟和带宽都比本地节点差得多：

本地节点访问延迟：~80ns
远程节点访问延迟：~140ns（1.75x）
本地带宽：~40GB/s per socket
远程带宽：~20GB/s per socket

并行排序的 NUMA 优化策略： 1. 使用 numactl 或 numa_alloc_onnode 确保数据分配在本地节点。 2. 分区排序时，确保每个线程处理本地节点上的数据。 3. 归并阶段尽量减少跨节点数据移动。

#include <numa.h>

/* NUMA 感知的数据分配 */
int *allocate_numa_aware(int n, int num_nodes) {
    int per_node = n / num_nodes;
    int *data = (int *)numa_alloc_interleaved(n * sizeof(int));

    /* 通过 first-touch 策略确保数据分布在正确的节点上 */
    #pragma omp parallel for schedule(static)
    for (int i = 0; i < n; i++) {
        data[i] = 0;  /* first touch 在对应线程的 NUMA 节点 */
    }
    return data;
}

十一、工程陷阱与最佳实践

11.1 常见陷阱表

11.2 选择并行排序算法的决策树

你的数据在哪里？
├─ 单机 CPU 内存
│   ├─ n < 100K → 串行排序（pdqsort）
│   ├─ 100K < n < 10M → 并行快排（TBB parallel_sort）
│   └─ n > 10M → 并行归并排序（考虑 NUMA）
├─ GPU 显存
│   ├─ 整数/浮点 → CUB DeviceRadixSort
│   └─ 自定义类型 → Thrust sort 或 ModernGPU merge sort
├─ 分布式集群
│   ├─ 均匀分布 → 采样排序
│   └─ 倾斜分布 → 自适应采样排序（增大 oversampling factor）
└─ 异构（CPU + GPU）
    └─ 数据量够大时 → GPU 排序 + CPU-GPU 流水线

11.3 关于稳定性

并行排序的稳定性是一个容易被忽略的问题。双调排序天然不稳定。并行快速排序也不稳定。只有并行归并排序可以保证稳定性，但前提是归并操作本身是稳定的（对相等元素保持来自左侧序列的在前）。

如果需要稳定排序又想用不稳定算法，标准技巧是把原始索引打包进 key：

struct StableKey {
    int value;
    int original_index;
    bool operator<(const StableKey& o) const {
        if (value != o.value) return value < o.value;
        return original_index < o.original_index;
    }
};

代价是 key 大小翻倍，对缓存不友好。

十二、个人观点与总结

12.1 我的一些看法

并行排序的 80/20 法则。在实际工程中，80% 的场景只需要两个方案：（1）多核 CPU 上用 std::sort 的并行版本（或 TBB），（2）GPU 上用 CUB 的基数排序。剩下 20% 的场景——自定义比较器、分布式环境、超大键值对——才需要深入理解本文讨论的各种算法。

排序网络的理论价值大于实用价值。双调排序的工程实现在 2024 年几乎不会出现在任何关键路径上——CUB 的基数排序在所有维度上都更好。但理解排序网络对理解并行计算的思维方式极有价值：如何将一个顺序依赖的问题拆解成固定的并行步骤，如何用 0-1 原理这样的技巧简化正确性证明。

GPU 排序的真正价值不在排序本身。在大多数 GPU 应用中（图形渲染、深度学习、物理模拟），排序是更大计算管道的一部分。数据已经在 GPU 显存中，排序只是中间步骤。这时 CPU-GPU 传输开销为零，GPU 排序的优势是碾压级的。那些”GPU 排序快 500 倍”的基准测试，在数据本来就在 GPU 上的场景中完全成立。

Merge Path 是我最欣赏的并行算法之一。它的优雅在于：用一个简单的二分查找，就将一个天然顺序的操作（归并）变成了完美并行的。没有复杂的数据结构，没有 tricky 的同步，只有一个精巧的数学观察。如果你只能记住本文的一个算法，记住 Merge Path。

不要迷信加速比。一个号称 100 倍加速的并行排序，如果单线程基准是未优化的冒泡排序，那这个数字毫无意义。始终用最优的单线程实现（pdqsort 级别）作为基准。

12.2 并行排序的未来

几个值得关注的方向：

1. CXL 和新内存架构。CXL（Compute Express Link）可能改变 NUMA 的拓扑，使得远程内存访问的延迟大幅降低。这会影响并行排序的数据放置策略。

2. 专用排序硬件。FPGA 上的排序网络可以跑到极高的频率。一些数据库加速卡已经集成了硬件排序单元。

3. 近内存计算。在内存颗粒旁边放计算单元，消除数据移动的瓶颈。这对排序这种内存密集型任务影响巨大。

4. 学习增强排序。Google 的 learned index 思路也被应用到排序——通过学习数据分布，更精确地分区。

12.3 关键公式汇总

附录：并行排序算法速查表

其中 w 是键的位宽，通常 32 或 64。

参考文献

1. K. E. Batcher, “Sorting Networks and Their Applications,” AFIPS Spring Joint Computer Conference, 1968.

2. M. Ajtai, J. Komlós, E. Szemerédi, “An O(n log n) Sorting Network,” STOC, 1983.

3. O. Green, R. McColl, D. Bader, “GPU Merge Path - A GPU Merging Algorithm,” ICS, 2012.

4. S. Baxter, “ModernGPU: Patterns and Idioms for GPU Computing,” 2013. https://moderngpu.github.io/

5. D. Merrill, “CUB: CUDA Unbound,” NVIDIA, 2015. https://nvlabs.github.io/cub/

6. H. Peters, O. Schulz-Hildebrandt, N. Luttenberger, “Fast In-Place, Comparison-Based Sorting with CUDA,” Parallel Computing, 2011.

7. A. Solomonik, L. Kale, “Highly Scalable Parallel Sorting,” IPDPS, 2010.

8. S. Sengupta, M. Harris, M. Garland, “Efficient Parallel Scan Algorithms for GPUs,” NVIDIA Technical Report, 2008.

9. J. Singler, P. Sanders, F. Putze, “MCSTL: The Multi-Core Standard Template Library,” Euro-Par, 2007.

10. Intel Threading Building Blocks Documentation. https://oneapi-src.github.io/oneTBB/

11. G. E. Blelloch, “Prefix Sums and Their Applications,” CMU Technical Report, 1990.

12. N. Satish, M. Harris, M. Garland, “Designing Efficient Sorting Algorithms for Manycore GPUs,” IPDPS, 2009.

算法系列导航：上一篇：外部排序 | 下一篇：排序基准测试

相关阅读：SIMD 算法设计模式 | pdqsort