10大规模并行处理器:高性能计算与CUDA编程(第4版)

并行模式:排序

GPU 上的排序

基数排序归并排序双调排序GPU排序优化

🎯 AI 时代的学习建议:不要只读理论。理解这一章的设计动机(为什么需要这个概念?)和边界条件(什么情况下会失效?)。当 AI 生成的代码涉及这些概念时,你能判断它是否用对了场景。AI 能帮你总结内容,但无法替代你亲手阅读、思考、建立直觉的过程。

第10章 并行模式:排序 - GPU上的排序 (Parallel Patterns: Sorting on the GPU)

一、导读

10.1 本章学习目标

通过本章的学习,读者将能够:

理解排序在并行计算中的重要性及其在GPU上实现的挑战

掌握Bitonic Sort算法的原理、实现和性能特征

深入理解Radix Sort(基数排序)的工作机制,包括直方图计算和前缀和应用

比较不同排序算法在GPU上的适用性,能够根据数据特点选择合适的算法

掌握共享内存在排序优化中的关键作用,减少全局内存访问

学会处理大规模数据的分块排序策略,协调多Block协作

了解现代GPU排序库(如Thrust、CUB)的实现特点和最佳实践

能够分析排序kernel的性能瓶颈,识别内存带宽和计算复杂度的影响

10.2 为什么本章如此重要

排序是计算机科学中最基础、最重要的操作之一。从数据库查询优化到图形渲染,从数据分析到机器学习,排序无处不在。在GPU上实现高效的排序具有特殊的意义和挑战。

首先,排序是一个计算密集型操作,传统的时间复杂度为O(N log N)。在GPU上,我们可以通过大规模并行化显著加速排序过程。然而,排序也涉及大量的数据移动和比较交换操作,这对内存子系统提出了很高的要求。

其次,不同的排序算法在GPU上的表现差异很大。经典的快速排序(QuickSort)由于其递归性和不规则的内存访问模式,在GPU上效率较低。而Bitonic Sort和Radix Sort等算法由于其规则的访问模式和高度并行的特性,更适合GPU架构。

第三,排序是许多高级算法的基础。基数排序内部使用直方图和前缀和;数据库的连接操作依赖排序来匹配记录;图形渲染中的深度排序需要高效的GPU排序。因此,掌握GPU排序技术对于构建高性能应用至关重要。

在实际项目中,低效的排序可能成为整个系统的瓶颈。例如,在粒子系统中,按深度排序粒子可以正确渲染透明物体;在数据库中,排序加速了范围查询和聚合操作;在机器学习中,排序用于特征选择和异常检测。因此,选择合适的高效排序算法对于性能优化至关重要。

本章将从经典的Bitonic Sort开始,逐步深入到更实用的Radix Sort,系统地介绍GPU排序的各种技术和优化策略。我们将通过详细的代码示例和性能分析,帮助读者建立起对GPU排序的全面理解。

10.3 核心问题引导

为什么传统的快速排序在GPU上表现不佳?GPU排序需要满足什么条件?

Bitonic Sort的工作原理是什么?它的并行复杂度是多少?

Radix Sort如何利用直方图和前缀和进行排序?它的时间复杂度是多少?

共享内存在排序中扮演什么角色?如何最大化数据复用?

如何处理超过单个Block处理能力的大规模数据排序?

Key-Value对的排序与纯Key排序有何不同?如何实现?

稳定排序和不稳定排序的区别是什么?哪些GPU排序算法是稳定的?

Thrust和CUB库中的排序实现有哪些优化技巧?何时应该使用库函数而非手写?


二、核心概念详解

10.1 GPU排序的挑战与要求

GPU排序的特殊要求:

在GPU上实现排序算法需要考虑以下因素:

规则的数据访问模式:GPU擅长处理规则的、可预测的访问模式。算法应尽量避免随机访问和指针追逐。

大规模并行性:算法应能充分利用GPU的数千个线程,避免串行瓶颈。

有限的分支分歧:排序涉及大量比较和条件交换,可能导致warp内的分支分歧。算法应尽量使相邻线程走相同的分支路径。

内存效率:排序需要频繁读写数据,内存访问模式直接影响性能。应优先使用共享内存,并确保全局内存的合并访问。

确定性执行时间:某些应用(如实时渲染)要求排序具有可预测的执行时间。算法的最坏情况复杂度应与平均情况相近。

为什么QuickSort不适合GPU:

快速排序是CPU上最常用的排序算法,平均时间复杂度为O(N log N)。但它在GPU上表现不佳,原因包括:

  • 递归结构:QuickSort的递归划分导致不规则的任务分解,难以映射到GPU的SIMT模型
  • 负载不均衡:分区后两个子数组的大小可能差异很大,导致线程利用率低
  • 不规则访问:分区操作涉及大量随机读写,缓存命中率低
  • 分支分歧:比较和交换操作导致warp内严重的分支分歧

相比之下,Bitonic Sort和Radix Sort具有规则的访问模式和确定的执行流程,更适合GPU。

10.2 Bitonic Sort算法

基本原理:

Bitonic Sort是一种基于比较网络的排序算法,由Ken Batcher于1968年提出。它的核心概念是bitonic序列:一个先单调递增后单调递减(或反之)的序列。

Bitonic Sort通过以下步骤工作:

构建Bitonic序列:将输入数组递归地划分为两半,一半升序排序,一半降序排序,合并形成bitonic序列。

Bitonic Merge:对于一个bitonic序列,通过比较相距N/2的元素对,将较小值放在前半部分,较大值放在后半部分。这一步产生两个较小的bitonic序列。

递归应用:对产生的两个较小bitonic序列递归应用bitonic merge,直到序列长度为1。

算法可视化(8元素示例):

Step 1: Create bitonic sequences of size 2
  [a0, a1] -> sort ascending/descending alternately

Step 2: Merge into bitonic sequences of size 4
  Compare elements at distance 2

Step 3: Merge into sorted sequence of size 8
  Compare elements at distance 4, then 2, then 1

并行复杂度:

  • 比较次数:O(N log²N)总比较次数
  • 并行步数:O(log²N)步
  • 处理器需求:每步需要N/2个比较器

Bitonic Sort的总工作量比QuickSort多一个log N因子,但其规则的访问模式和完全并行的特性使其在GPU上具有竞争力。

CUDA实现:

cpp
__global__ void bitonicSort(float* data, int N, int stage, int pass) {
    int tid = threadIdx.x + blockIdx.x * blockDim.x;
    int pairDistance = 1 << (stage - pass);
    int blockWidth = 2 * pairDistance;
    
    int leftId = (tid / pairDistance) % 2 == 0 
                 ? tid - (tid % pairDistance) 
                 : tid + pairDistance;
    
    if (leftId >= N || leftId < 0) return;
    
    // Determine sort direction based on block index
    int blockIdx_local = tid / blockWidth;
    bool ascending = (blockIdx_local % 2) == 0;
    
    float leftVal = data[leftId];
    float rightVal = data[leftId + pairDistance];
    
    if ((ascending && leftVal > rightVal) || 
        (!ascending && leftVal < rightVal)) {
        data[leftId] = rightVal;
        data[leftId + pairDistance] = leftVal;
    }
}

// Host function to launch all stages
void gpuBitonicSort(float* d_data, int N) {
    // N must be power of 2
    int threadsPerBlock = N / 2;
    if (threadsPerBlock > 512) threadsPerBlock = 512;
    
    for (int stage = 1; stage < (int)(log2((double)N) + 1); stage++) {
        for (int pass = 0; pass < stage; pass++) {
            bitonicSort<<<(N + threadsPerBlock * 2 - 1) / (threadsPerBlock * 2), 
                          threadsPerBlock>>>(d_data, N, stage, pass);
            cudaDeviceSynchronize();
        }
    }
}

这个实现的缺点是每次迭代都需要单独的kernel启动和同步,开销较大。优化的版本可以在单个kernel内完成所有阶段,使用__syncthreads()进行Block内同步。

优化:单Kernel Bitonic Sort:

cpp
__global__ void bitonicSortSingleKernel(float* data, int N) {
    extern __shared__ float shared[];
    
    int tid = threadIdx.x;
    int globalIdx = blockIdx.x * blockDim.x + threadIdx.x;
    
    // Load data into shared memory
    if (globalIdx < N) {
        shared[tid] = data[globalIdx];
    }
    __syncthreads();
    
    int blockSize = blockDim.x;
    
    // Bitonic sort within shared memory
    for (int stage = 1; stage < blockSize; stage *= 2) {
        int pairDistance = stage;
        int blockWidth = 2 * stage;
        
        int leftId = (tid / pairDistance) % 2 == 0 
                     ? (tid / blockWidth) * blockWidth + (tid % pairDistance)
                     : (tid / blockWidth) * blockWidth + pairDistance + (tid % pairDistance);
        
        int rightId = leftId ^ pairDistance;  // XOR flips the relevant bit
        
        if (rightId >= blockSize) continue;
        
        bool ascending = ((tid / blockWidth) % 2) == 0;
        
        float leftVal = shared[leftId];
        float rightVal = shared[rightId];
        
        if ((ascending && leftVal > rightVal) || 
            (!ascending && leftVal < rightVal)) {
            shared[leftId] = rightVal;
            shared[rightId] = leftVal;
        }
        __syncthreads();
    }
    
    // Write back
    if (globalIdx < N) {
        data[globalIdx] = shared[tid];
    }
}

这个版本将数据加载到共享内存中,在Block内进行完整的bitonic sort,减少了全局内存访问。但它仅限于Block大小的数据。

10.3 Radix Sort算法

基本原理:

Radix Sort是一种非比较排序算法,它通过逐位(digit-by-digit)处理来排序数据。对于整数或浮点数,我们可以将其二进制表示划分为多个位段(如每次处理4位或8位),然后对每个位段进行稳定排序。

Radix Sort的核心步骤:

提取位段:对于当前处理的位段,从每个key中提取对应的位值。

直方图计算:统计每个位值(0到基数-1)的出现次数。

前缀和(Scan):对直方图进行exclusive scan,得到每个位值的起始位置。

分散(Scatter):根据scan结果,将元素放置到输出数组的正确位置。

重复:对下一个位段重复上述步骤,直到所有位段处理完毕。

基数选择:

基数(radix)决定了每次处理的位数。常见的选择包括:

  • Radix 2⁴ = 16:每次处理4位,需要8次迭代(对于32位整数)
  • Radix 2⁸ = 256:每次处理8位,需要4次迭代
  • Radix 2¹⁶ = 65536:每次处理16位,需要2次迭代

较大的基数减少了迭代次数,但增加了直方图的大小和scan的开销。通常,radix 256是一个较好的平衡点。

复杂度分析:

  • 时间复杂度:O(k × N),其中k是位数除以每次处理的位数
  • 空间复杂度:O(N + radix),需要额外的输出数组和直方图
  • 稳定性:Radix Sort是稳定排序,保持相等元素的原始顺序

对于32位整数,使用radix 256时,k = 4,时间复杂度为O(4N) = O(N),优于比较排序的O(N log N)。

CUDA实现:

cpp
#define RADIX_BITS 8
#define RADIX_SIZE (1 << RADIX_BITS)  // 256
#define RADIX_MASK (RADIX_SIZE - 1)

// Kernel to compute histogram per block
__global__ void radixHistogram(const unsigned int* keys, 
                                unsigned int* histograms,
                                int N, int bitOffset) {
    extern __shared__ unsigned int sharedHist[RADIX_SIZE];
    
    int tid = threadIdx.x;
    int globalIdx = blockIdx.x * blockDim.x + threadIdx.x;
    
    // Initialize shared histogram
    if (tid < RADIX_SIZE) {
        sharedHist[tid] = 0;
    }
    __syncthreads();
    
    // Each thread processes multiple elements
    for (int i = globalIdx; i < N; i += blockDim.x * gridDim.x) {
        unsigned int key = keys[i];
        unsigned int digit = (key >> bitOffset) & RADIX_MASK;
        atomicAdd(&sharedHist[digit], 1);
    }
    __syncthreads();
    
    // Write shared histogram to global histogram
    // Each thread writes multiple bins
    for (int i = tid; i < RADIX_SIZE; i += blockDim.x) {
        histograms[blockIdx.x * RADIX_SIZE + i] = sharedHist[i];
    }
}

// Kernel to perform radix sort step (histogram + scan + scatter)
__global__ void radixSortStep(const unsigned int* inputKeys,
                               unsigned int* outputKeys,
                               const unsigned int* values,  // Optional values array
                               unsigned int* outputValues,
                               unsigned int* offsets,       // Pre-computed scan results
                               int N, int bitOffset) {
    int tid = threadIdx.x;
    int globalIdx = blockIdx.x * blockDim.x + threadIdx.x;
    
    // Use shared memory for local histogram and prefix sums
    extern __shared__ unsigned int shared[];
    unsigned int* localHist = shared;
    unsigned int* localScan = shared + RADIX_SIZE;
    
    // Initialize
    if (tid < RADIX_SIZE) {
        localHist[tid] = 0;
    }
    __syncthreads();
    
    // Compute local histogram
    for (int i = globalIdx; i < N; i += blockDim.x * gridDim.x) {
        unsigned int digit = (inputKeys[i] >> bitOffset) & RADIX_MASK;
        atomicAdd(&localHist[digit], 1);
    }
    __syncthreads();
    
    // Compute exclusive scan of local histogram
    // (Using simple sequential scan for clarity)
    if (tid == 0) {
        localScan[0] = 0;
        for (int i = 1; i < RADIX_SIZE; i++) {
            localScan[i] = localScan[i - 1] + localHist[i - 1];
        }
    }
    __syncthreads();
    
    // Add global offsets
    if (tid < RADIX_SIZE) {
        localScan[tid] += offsets[blockIdx.x * RADIX_SIZE + tid];
    }
    __syncthreads();
    
    // Scatter elements to output
    for (int i = globalIdx; i < N; i += blockDim.x * gridDim.x) {
        unsigned int key = inputKeys[i];
        unsigned int digit = (key >> bitOffset) & RADIX_MASK;
        
        // Atomic add to get unique position within this digit bin
        unsigned int pos = atomicAdd(&localScan[digit], 1);
        
        outputKeys[pos] = key;
        if (values != NULL) {
            outputValues[pos] = values[i];
        }
    }
}

// Host function for full radix sort
void gpuRadixSort(unsigned int* d_keys, unsigned int* d_values, int N) {
    unsigned int *d_tempKeys, *d_tempValues;
    cudaMalloc((void**)&d_tempKeys, N * sizeof(unsigned int));
    cudaMalloc((void**)&d_tempValues, N * sizeof(unsigned int));
    
    unsigned int *currentKeys = d_keys;
    unsigned int *currentValues = d_values;
    unsigned int *nextKeys = d_tempKeys;
    unsigned int *nextValues = d_tempValues;
    
    int blockSize = 256;
    int gridSize = (N + blockSize - 1) / blockSize;
    
    // Allocate offsets array
    unsigned int *d_offsets;
    cudaMalloc((void**)&d_offsets, gridSize * RADIX_SIZE * sizeof(unsigned int));
    
    // Process each bit segment
    for (int bitOffset = 0; bitOffset < 32; bitOffset += RADIX_BITS) {
        // Step 1: Compute global histogram
        // (Omitted for brevity - use CUB or thrust for production code)
        
        // Step 2: Exclusive scan of histogram
        // (Omitted for brevity)
        
        // Step 3: Scatter
        radixSortStep<<<gridSize, blockSize, 
                         (RADIX_SIZE * 2) * sizeof(unsigned int)>>>(
            currentKeys, nextKeys, currentValues, nextValues,
            d_offsets, N, bitOffset);
        
        // Swap buffers
        unsigned int* temp = currentKeys;
        currentKeys = nextKeys;
        nextKeys = temp;
        
        temp = currentValues;
        currentValues = nextValues;
        nextValues = temp;
    }
    
    // If final result is in temporary buffer, copy back
    if (currentKeys != d_keys) {
        cudaMemcpy(d_keys, currentKeys, N * sizeof(unsigned int), 
                   cudaMemcpyDeviceToDevice);
        if (d_values != NULL) {
            cudaMemcpy(d_values, currentValues, N * sizeof(unsigned int), 
                       cudaMemcpyDeviceToDevice);
        }
    }
    
    cudaFree(d_tempKeys);
    cudaFree(d_tempValues);
    cudaFree(d_offsets);
}

这个实现展示了Radix Sort的基本框架。在实际应用中,直方图计算和scan可以使用CUB库的高效实现。

10.4 共享内存在排序中的优化

减少全局内存访问:

排序涉及大量的数据读写,将数据暂存到共享内存可以显著减少全局内存访问。在Bitonic Sort中,我们已经在Block内使用共享内存。在Radix Sort中,共享内存用于:

局部直方图:每个Block维护自己的直方图,最后合并到全局直方图

局部Scan:在共享内存中进行scan,减少全局同步

数据缓冲:暂存待排序的数据,批量写入全局内存

Bank Conflict避免:

共享内存的bank conflict在排序中尤其常见,因为多个线程可能同时更新同一个bin的计数器。使用原子操作可以正确处理竞争,但会降低性能。优化方法包括:

私有直方图:每个线程维护私有直方图,最后归约

Warp级归约:在warp内先归约计数,再原子更新共享直方图

填充:在直方图数组中添加填充,避免多个bin映射到同一bank

10.5 大规模数据的分块排序

当数据规模超过单个Block的处理能力时,需要进行分块排序。

策略1:多阶段Radix Sort

将数据划分为多个Block大小的段

每个Block独立计算局部直方图

合并局部直方图,计算全局前缀和

根据全局前缀和,将元素分散到正确的位置

策略2:混合排序

对于非常大的数据集,可以结合多种排序算法:

使用Radix Sort对数据进行粗粒度排序(如只处理高16位)

对每个桶使用Bitonic Sort或其他方法进行细粒度排序

策略3:外部排序

当数据超出显存容量时,需要使用外部排序:

将数据分批加载到GPU

对每批数据进行GPU排序

将排序后的批次写回主机内存

使用多路归并合并所有批次

10.6 Key-Value排序

在许多应用中,我们需要根据Key对Key-Value对进行排序,同时保持Key和Value的对应关系。

实现方法:

在Radix Sort的scatter阶段,同时移动Key和Value:

cpp
outputKeys[pos] = key;
outputValues[pos] = values[i];

这确保了Key和Value始终一起移动,保持对应关系。

优化:

如果Value较大(如结构体),移动Value的开销可能很大。可以考虑:

索引排序:只对索引排序,然后通过索引间接访问Value

分离存储:将Key和Value存储在不同的数组中,只移动必要的部分

10.7 稳定排序 vs 不稳定排序

稳定排序:相等元素的相对顺序保持不变。Radix Sort是稳定排序。

不稳定排序:相等元素的相对顺序可能改变。Bitonic Sort是不稳定排序。

在某些应用中(如多关键字排序),稳定性很重要。Radix Sort天然支持稳定排序,而Bitonic Sort需要额外处理才能保持稳定。

10.8 使用Thrust和CUB库

在生产代码中,推荐使用经过充分优化的库函数:

Thrust:

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

thrust::device_vector<int> d_data(N);
// ... fill data ...
thrust::sort(d_data.begin(), d_data.end());

CUB:

cpp
#include <cub/cub.cuh>

void* d_tempStorage = NULL;
size_t tempStorageBytes = 0;

cub::DeviceRadixSort::SortKeys(d_tempStorage, tempStorageBytes, 
                                d_keys, d_sortedKeys, N);
cudaMalloc(&d_tempStorage, tempStorageBytes);
cub::DeviceRadixSort::SortKeys(d_tempStorage, tempStorageBytes, 
                                d_keys, d_sortedKeys, N);
cudaFree(d_tempStorage);

这些库提供了高度优化的实现,支持各种数据类型和排序选项。


三、程序员的视角

3.1 在实际项目中的应用

数据库:排序用于ORDER BY查询、索引构建、连接操作等。GPU加速的数据库引擎使用Radix Sort加速这些操作。

图形渲染:在粒子系统渲染中,按深度排序粒子可以正确渲染透明物体。在光线追踪中,排序射线可以提高缓存命中率。

数据分析:排序用于排名、百分位计算、去重等操作。在大数据处理中,GPU排序可以显著加速这些任务。

机器学习:在某些特征选择算法中,需要按重要性排序特征。K近邻算法也需要排序来找到最近的邻居。

3.2 常见陷阱与误区

数据大小不是基数的倍数:Radix Sort要求数据大小适合基数划分。对于剩余部分需要特殊处理。

忘记交换缓冲区:Radix Sort的每次迭代需要在输入和输出缓冲区之间交换。忘记交换会导致错误结果。

溢出问题:直方图计数可能溢出32位整数,特别是对于大规模数据。使用64位整数或检查溢出。

稳定性假设:Bitonic Sort是不稳定的,如果需要稳定排序,应选择Radix Sort。

忽视Key-Value关联:排序Key时忘记同步移动Value,导致数据错位。

3.3 与其他层级的联系

与前缀和的联系:Radix Sort的核心是直方图的前缀和。高效的scan实现直接提升排序性能。

与归约的联系:直方图计算本质上是一种归约。归约优化技术(如共享内存归约)也适用于直方图。

与内存层次的联系:排序是内存密集型操作,共享内存和缓存的利用对性能至关重要。


四、代码实战

4.1 完整的Radix Sort实现

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

#define RADIX_BITS 8
#define RADIX_SIZE (1 << RADIX_BITS)
#define BLOCK_SIZE 256

// Simplified radix sort for demonstration
__global__ void radixSortPass(const unsigned int* input, unsigned int* output,
                               int N, int bitOffset) {
    int idx = blockIdx.x * blockDim.x + threadIdx.x;
    
    if (idx < N) {
        unsigned int key = input[idx];
        unsigned int digit = (key >> bitOffset) & (RADIX_SIZE - 1);
        
        // In practice, use shared memory histogram and scan
        // This is a simplified version
        output[idx] = key;  // Placeholder
    }
}

int main() {
    int N = 1000000;
    unsigned int *h_data = (unsigned int*)malloc(N * sizeof(unsigned int));
    
    for (int i = 0; i < N; i++) {
        h_data[i] = rand();
    }
    
    unsigned int *d_data, *d_temp;
    cudaMalloc((void**)&d_data, N * sizeof(unsigned int));
    cudaMalloc((void**)&d_temp, N * sizeof(unsigned int));
    
    cudaMemcpy(d_data, h_data, N * sizeof(unsigned int), 
               cudaMemcpyHostToDevice);
    
    printf("Radix Sort on GPU demonstrated.\n");
    printf("For production use, consider CUB or Thrust libraries.\n");
    
    cudaFree(d_data);
    cudaFree(d_temp);
    free(h_data);
    
    return 0;
}

4.2 使用Thrust进行排序

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

int main() {
    const int N = 1000000;
    
    thrust::device_vector<int> d_data(N);
    
    // Fill with random data
    for (int i = 0; i < N; i++) {
        d_data[i] = rand();
    }
    
    // Sort
    thrust::sort(d_data.begin(), d_data.end());
    
    // Verify
    bool sorted = true;
    for (int i = 1; i < N; i++) {
        if (d_data[i] < d_data[i-1]) {
            sorted = false;
            break;
        }
    }
    
    printf("Thrust sort: %s\n", sorted ? "PASSED" : "FAILED");
    
    return 0;
}

五、深入思考

5.1 历史背景

Bitonic Sort由Ken Batcher于1968年提出,是最早的并行排序网络之一。Radix Sort的历史更早,可以追溯到打孔卡时代。

在GPU领域,早期的研究集中在Bitonic Sort上,因为其规则的访问模式。随着GPU计算能力的提升,Radix Sort因其更好的时间复杂度而受到青睐。NVIDIA的CUB库和Thrust库提供了生产质量的实现。

5.2 前沿方向

Sample Sort:一种自适应的并行排序算法,根据数据分布动态调整分区。

GPU-Aware排序:利用新一代GPU的硬件特性(如异步拷贝、Tensor Core)加速排序。

外部排序优化:针对超出显存容量的超大数据集,研究高效的多路归并策略。

5.3 延伸阅读

NVIDIA CUB Library Documentation - DeviceRadixSort

Thrust Library Documentation - sort

Satish et al., "Designing Efficient Sorting Algorithms for Manycore GPUs", IPDPS 2009

Batcher, "Sorting Networks and Their Applications", 1968


六、实验笔记

实验1:Bitonic Sort vs Radix Sort性能对比

目的:比较两种排序算法在不同数据规模下的性能。

步骤:

实现两种算法的CUDA版本

测试N = 10^4, 10^5, ..., 10^7的执行时间

分析时间复杂度差异

预期:对于小数据,Bitonic Sort可能更快,因为常数因子小。对于大数据,Radix Sort的O(N)复杂度优势显现。

实验2:共享内存优化效果

目的:量化共享内存对Radix Sort直方图计算的加速效果。

步骤:

实现直接使用全局内存原子操作的版本

实现使用共享内存局部直方图的版本

比较两者的执行时间和内存事务数量

预期:共享内存版本应显著减少全局内存原子操作,提高带宽利用率。


七、本章小结

排序是GPU计算中的重要操作,本章介绍了两种主要的GPU排序算法:Bitonic Sort和Radix Sort。

Bitonic Sort是一种基于比较网络的排序算法,具有规则的访问模式和确定的执行时间。它的总工作量为O(N log²N),适合中小规模数据的Block内排序。

Radix Sort是一种非比较排序算法,通过逐位处理和稳定散射实现排序。它的时间复杂度为O(kN),其中k是位数除以每次处理的位数。对于32位整数,使用radix 256时,k=4,实际性能优异。Radix Sort是稳定排序,支持Key-Value对排序。

共享内存在排序中扮演关键角色,用于局部直方图、scan和数据缓冲,显著减少全局内存访问。对于大规模数据,需要采用分块策略,协调多Block协作。

在生产代码中,推荐使用Thrust或CUB库提供的排序实现,它们经过充分优化,可以达到接近硬件峰值的性能。


八、习题

Bitonic Sort追踪:手动追踪Bitonic Sort对8元素数组[5, 2, 8, 1, 9, 3, 7, 4]的排序过程。

Radix Sort迭代次数:对于64位整数,使用radix 256,需要多少次迭代?使用radix 16呢?

稳定性分析:为什么Radix Sort是稳定排序?Bitonic Sort为什么不稳定?

直方图溢出:如果N = 10^9,使用32位整数存储直方图bin,是否可能溢出?如何避免?

Key-Value排序:修改Radix Sort kernel,支持同时对Key和Value排序。

反向排序:如何修改Radix Sort以实现降序排序?

浮点数排序:IEEE 754浮点数的二进制表示不能直接用于Radix Sort。如何预处理浮点数使其可以正确排序?

多关键字排序:设计一个算法,先按Key1排序,再按Key2稳定排序。说明如何使用Radix Sort实现。

性能建模:估算Radix Sort的全局内存访问次数。与QuickSort相比,哪个的内存效率更高?

综合实现:实现一个完整的Radix Sort,包括直方图计算、scan和scatter。使用CUB库优化scan部分。测量排序10^7个整数的时间。

10/15并行模式
0%