第8章并行模式：归约 - 从数据到结果 (Parallel Patterns: Reduction - From Data to Result)

一、导读

8.1 本章学习目标

通过本章的学习，读者将能够：

深入理解归约（Reduction）操作的定义、特性和在并行计算中的核心地位

掌握树形归约的基本原理，能够设计高效的并行归约算法

精通共享内存在归约优化中的应用，包括交错寻址和银行冲突避免

理解并实现多种归约优化技术：展开最后warp、多元素处理、原子操作替代

掌握并行前缀和与归约的关系，了解扫描操作的基础

学会处理大规模数据的分段归约策略，协调Block级和Grid级归约

能够分析归约kernel的性能瓶颈，识别内存带宽和同步开销的影响

了解归约在其他并行算法中的应用，如直方图、argmax/argmin等变体

8.2 为什么本章如此重要

归约是并行计算中最基础、最重要的模式之一。它将一个数组通过二元运算符（如加法、乘法、最大值、最小值等）缩减为单个值。看似简单的操作，在并行环境下却有着丰富的优化空间。

归约的重要性体现在多个方面：

首先，归约是许多复杂算法的核心子程序。求和用于统计分析和数值积分；最大值/最小值用于数据范围确定和异常检测；逻辑与/或用于条件聚合；甚至矩阵乘法和卷积内部也包含归约操作。可以说，掌握了归约，就掌握了并行编程的半壁江山。

其次，归约是学习并行优化的绝佳案例。它涉及共享内存使用、线程同步、树形并行、warp级优化、原子操作等多个关键技术点。通过对归约的逐步优化，我们可以清晰地看到每个优化技术的贡献，建立起系统的性能调优思维。

第三，归约的性能特征非常典型：它是内存带宽受限的操作，因为每个数据元素只参与一次运算，计算强度很低。这使得归约成为测试和优化内存子系统的理想基准。同时，归约需要同步机制来协调部分结果的合并，这引入了同步开销的考量。

在实际项目中，低效的归约实现可能成为整个应用的瓶颈。例如，在深度学习训练中，梯度归约需要在多GPU之间同步；在科学模拟中，能量和动量的守恒检查需要全局归约；在数据分析中，聚合查询依赖高效的归约。因此，掌握归约优化技术对于构建高性能应用至关重要。

本章将从最简单的串行归约开始，逐步引入并行化策略和各种优化技术，最终达到接近硬件理论峰值的性能。我们将通过详细的代码示例和性能测量，展示每一步优化的效果和原理。

8.3 核心问题引导

什么是归约操作？它需要满足什么数学性质才能正确并行化？

树形归约的工作原理是什么？为什么它的并行复杂度是O(log N)？

如何使用共享内存实现Block级归约？交错寻址和相邻寻址有什么区别？

什么是"展开最后warp"优化？它为什么能消除同步开销？

如何让每个线程处理多个元素以提高效率？这对occupancy有何影响？

当数据规模超过一个Block的处理能力时，如何进行多Block归约？

原子操作能否替代归约？在什么情况下原子操作更高效？

归约的计算强度是多少？为什么它通常是内存带宽受限的？

二、核心概念详解

8.1 归约操作的定义与性质

归约的定义：

归约（Reduction）是将一个数据集合通过二元运算符⊕合并为单个值的操作。数学上表示为：

$$result = x_0 \oplus x_1 \oplus x_2 \oplus ... \oplus x_{N-1}$$

其中⊕是二元结合运算符。

常见的归约运算符：

运算符	初始值	应用
加法 (+)	0	求和、平均值
乘法 (×)	1	乘积、阶乘
最大值 (max)	-∞	找最大值、argmax
最小值 (min)	+∞	找最小值、argmin
逻辑与 (&&)	true	全真判断
逻辑或 (\	\	)	false	存在性判断
按位与 (&)	全1	掩码合并
按位或 (\	)	0	标志合并

并行化的数学要求：

为了正确并行化归约，运算符必须满足：

结合律（Associativity）：$(a \oplus b) \oplus c = a \oplus (b \oplus c)$

结合律保证了我们可以任意分组数据进行并行计算，而不影响最终结果。加法、乘法、max、min都满足结合律。

交换律（Commutativity）（可选但推荐）：$a \oplus b = b \oplus a$

交换律不是严格必需的，但它允许我们任意重排数据顺序。浮点加法在理论上不满足交换律（由于舍入误差），但在实践中通常当作可交换处理。

注意：浮点运算由于舍入误差，不同的归约顺序可能产生略微不同的结果。这在数值敏感的应用中需要考虑。

8.2 树形归约算法

基本思想：

树形归约通过将数据组织成二叉树结构，自底向上逐层合并。每一层将相邻的两个元素合并为一个，直到只剩下一个元素。

Level 0: [1] [2] [3] [4] [5] [6] [7] [8]
          \ /   \ /   \ /   \ /
Level 1:  [3]   [7]   [11]  [13]
           \     /     \     /
Level 2:    [10]       [24]
              \         /
Level 3:        [34]

对于一个大小为N的数组，树形归约需要log₂N层，每层的并行工作量减半。总的工作量仍然是N-1次运算，但并行深度从串行的N-1降为log₂N。

并行复杂度：

时间复杂度：O(log N)并行步数
工作复杂度：O(N)总运算次数（与串行相同）
处理器需求：最多需要N/2个处理器（在第一层）

在实际GPU实现中，我们受限于Block内的线程数量，通常使用分层归约：先在Block内归约，然后在全局级别合并各Block的结果。

8.3 基于共享内存的Block级归约

基本实现：

使用共享内存存储Block内线程的部分结果，然后通过树形结构逐步归约。

cpp

__global__ void reductionBasic(const float* input, float* output, int N) {
    extern __shared__ float shared[];
    
    int tid = threadIdx.x;
    int idx = blockIdx.x * blockDim.x + threadIdx.x;
    
    // Load input into shared memory
    shared[tid] = (idx < N) ? input[idx] : 0.0f;
    __syncthreads();
    
    // Tree-based reduction
    for (int stride = blockDim.x / 2; stride > 0; stride >>= 1) {
        if (tid < stride) {
            shared[tid] += shared[tid + stride];
        }
        __syncthreads();
    }
    
    // First thread writes block result
    if (tid == 0) {
        atomicAdd(output, shared[0]);
    }
}

这个实现的问题是每个步长都需要一次__syncthreads()，同步开销较大。

交错寻址 vs 相邻寻址：

上面的实现使用交错寻址（interleaved addressing），即thread i合并shared[i]和shared[i + stride]。这种方法在早期步长中会导致bank conflict，因为访问的地址间隔较大。

相邻寻址（sequential addressing）让相邻线程处理相邻数据：

cpp

__global__ void reductionAdjacent(const float* input, float* output, int N) {
    extern __shared__ float shared[];
    
    int tid = threadIdx.x;
    int idx = blockIdx.x * (blockDim.x * 2) + threadIdx.x;
    
    // Each thread loads two elements
    float val = (idx < N) ? input[idx] : 0.0f;
    if (idx + blockDim.x < N) {
        val += input[idx + blockDim.x];
    }
    shared[tid] = val;
    __syncthreads();
    
    // Reduction with sequential addressing
    for (int stride = blockDim.x / 2; stride > 0; stride >>= 1) {
        if (tid < stride) {
            shared[tid] += shared[tid + stride];
        }
        __syncthreads();
    }
    
    if (tid == 0) {
        atomicAdd(output, shared[0]);
    }
}

相邻寻址在早期步长中具有更好的内存访问模式，因为相邻线程访问相邻地址。

8.4 优化技术详解

优化1：展开最后warp

在归约的最后阶段，当stride ≤ 32（warp大小）时，所有活跃线程都在同一个warp内。由于warp内的线程天然同步（在Volta之前），可以省略__syncthreads()，手动展开循环以减少分支开销。

cpp

__device__ float warpReduce(volatile float* shared, int tid) {
    shared[tid] += shared[tid + 32];
    shared[tid] += shared[tid + 16];
    shared[tid] += shared[tid + 8];
    shared[tid] += shared[tid + 4];
    shared[tid] += shared[tid + 2];
    shared[tid] += shared[tid + 1];
    return shared[tid];
}

__global__ void reductionUnrolledWarp(const float* input, float* output, int N) {
    extern __shared__ float shared[];
    
    int tid = threadIdx.x;
    int idx = blockIdx.x * (blockDim.x * 2) + threadIdx.x;
    
    float val = (idx < N) ? input[idx] : 0.0f;
    if (idx + blockDim.x < N) {
        val += input[idx + blockDim.x];
    }
    shared[tid] = val;
    __syncthreads();
    
    // Standard reduction until stride <= 32
    for (int stride = blockDim.x / 2; stride > 32; stride >>= 1) {
        if (tid < stride) {
            shared[tid] += shared[tid + stride];
        }
        __syncthreads();
    }
    
    // Unroll last warp
    if (tid < 32) {
        volatile float* vshared = shared;  // volatile prevents compiler optimization
        warpReduce(vshared, tid);
    }
    
    if (tid == 0) {
        atomicAdd(output, shared[0]);
    }
}

这个优化消除了最后5次同步，显著减少了同步开销。使用volatile关键字防止编译器优化掉对共享内存的重复访问。

优化2：完全展开归约树

对于固定的Block大小（如256），可以完全展开归约循环，避免循环控制开销：

cpp

template <unsigned int blockSize>
__device__ void reduceStep(volatile float* shared, int tid) {
    if (blockSize >= 64) shared[tid] += shared[tid + 32];
    if (blockSize >= 32) shared[tid] += shared[tid + 16];
    if (blockSize >= 16) shared[tid] += shared[tid + 8];
    if (blockSize >= 8)  shared[tid] += shared[tid + 4];
    if (blockSize >= 4)  shared[tid] += shared[tid + 2];
    if (blockSize >= 2)  shared[tid] += shared[tid + 1];
}

template <unsigned int blockSize>
__global__ void reductionFullyUnrolled(const float* input, float* output, int N) {
    extern __shared__ float shared[];
    
    int tid = threadIdx.x;
    int idx = blockIdx.x * (blockSize * 2) + threadIdx.x;
    
    float val = (idx < N) ? input[idx] : 0.0f;
    if (idx + blockSize < N) {
        val += input[idx + blockSize];
    }
    shared[tid] = val;
    __syncthreads();
    
    if (blockSize >= 512 && tid < 256) {
        shared[tid] += shared[tid + 256];
    }
    __syncthreads();
    if (blockSize >= 256 && tid < 128) {
        shared[tid] += shared[tid + 128];
    }
    __syncthreads();
    if (blockSize >= 128 && tid < 64) {
        shared[tid] += shared[tid + 64];
    }
    __syncthreads();
    
    if (tid < 32) {
        volatile float* vshared = shared;
        reduceStep<blockSize>(vshared, tid);
    }
    
    if (tid == 0) {
        atomicAdd(output, shared[0]);
    }
}

通过模板参数指定Block大小，编译器可以在编译时确定哪些条件为真，消除不必要的分支。

优化3：多元素处理

让每个线程处理多个元素，可以减少Block数量，从而减少最后的原子操作竞争。

cpp

__global__ void reductionMultiElement(const float* input, float* output, int N) {
    extern __shared__ float shared[];
    
    int tid = threadIdx.x;
    int gridSize = blockDim.x * gridDim.x;
    int idx = blockIdx.x * blockDim.x + threadIdx.x;
    
    // Each thread processes multiple elements across the entire array
    float val = 0.0f;
    for (int i = idx; i < N; i += gridSize) {
        val += input[i];
    }
    shared[tid] = val;
    __syncthreads();
    
    // Standard reduction
    for (int stride = blockDim.x / 2; stride > 32; stride >>= 1) {
        if (tid < stride) {
            shared[tid] += shared[tid + stride];
        }
        __syncthreads();
    }
    
    if (tid < 32) {
        volatile float* vshared = shared;
        if (blockDim.x >= 64) vshared[tid] += vshared[tid + 32];
        if (blockDim.x >= 32) vshared[tid] += vshared[tid + 16];
        if (blockDim.x >= 16) vshared[tid] += vshared[tid + 8];
        if (blockDim.x >= 8)  vshared[tid] += vshared[tid + 4];
        if (blockDim.x >= 4)  vshared[tid] += vshared[tid + 2];
        if (blockDim.x >= 2)  vshared[tid] += vshared[tid + 1];
    }
    
    if (tid == 0) {
        atomicAdd(output, shared[0]);
    }
}

这种方法的优点是只需要一个kernel启动，不需要分段的Grid级归约。缺点是当N很大时，每个线程需要处理很多元素，可能导致负载不均衡。

优化4：使用原子操作简化

对于某些归约（特别是加法），可以直接使用原子操作，避免复杂的归约树：

cpp

__global__ void reductionAtomic(const float* input, float* output, int N) {
    int idx = blockIdx.x * blockDim.x + threadIdx.x;
    if (idx < N) {
        atomicAdd(output, input[idx]);
    }
}

这种方法极其简单，但性能取决于原子操作的竞争程度。当大量线程同时原子写入同一地址时，序列化会导致性能下降。但对于小规模数据或作为多级归约的最后一步，原子操作是一个实用的选择。

在Ampere及以后的架构中，atomicAdd对FP32的支持得到了硬件加速，性能显著提升。

8.5 多Block归约策略

当数据规模超过单个Block的处理能力时，需要进行多Block归约。常见策略有：

策略1：两阶段归约

第一阶段：每个Block计算局部归约结果，写入全局数组

第二阶段：对局部结果数组进行第二次归约

cpp

// Phase 1: Block-level reduction
__global__ void reductionPhase1(const float* input, float* partialSums, int N) {
    extern __shared__ float shared[];
    
    int tid = threadIdx.x;
    int idx = blockIdx.x * blockDim.x + threadIdx.x;
    
    float val = 0.0f;
    for (int i = idx; i < N; i += blockDim.x * gridDim.x) {
        val += input[i];
    }
    shared[tid] = val;
    __syncthreads();
    
    // ... reduction tree ...
    
    if (tid == 0) {
        partialSums[blockIdx.x] = shared[0];
    }
}

// Phase 2: Reduce partial sums (can use same kernel or CPU)
void multiBlockReduction(const float* input, float* output, int N) {
    int blockSize = 256;
    int gridSize = min(1024, (N + blockSize * 2 - 1) / (blockSize * 2));
    
    float *d_partialSums;
    cudaMalloc((void**)&d_partialSums, gridSize * sizeof(float));
    
    // Phase 1
    reductionPhase1<<<gridSize, blockSize, blockSize * sizeof(float)>>>(
        input, d_partialSums, N);
    
    // Phase 2: reduce partial sums on GPU or CPU
    float *h_partialSums = (float*)malloc(gridSize * sizeof(float));
    cudaMemcpy(h_partialSums, d_partialSums, gridSize * sizeof(float), 
               cudaMemcpyDeviceToHost);
    
    float sum = 0.0f;
    for (int i = 0; i < gridSize; i++) {
        sum += h_partialSums[i];
    }
    *output = sum;
    
    cudaFree(d_partialSums);
    free(h_partialSums);
}

策略2：单Kernel多Pass

使用CUDA动态并行或全局原子操作，在一个Kernel内完成所有归约。这种方法减少了kernel启动开销，但实现更复杂。

策略3：使用CUB库

NVIDIA的CUB库提供了高度优化的归约实现，支持各种数据类型和运算符。在生产代码中，推荐使用CUB而非手写归约。

cpp

#include <cub/cub.cuh>

void cubReduction(const float* d_input, float* d_output, int N) {
    void* d_tempStorage = NULL;
    size_t tempStorageBytes = 0;
    
    // Determine temporary device storage requirements
    cub::DeviceReduce::Sum(d_tempStorage, tempStorageBytes, 
                           d_input, d_output, N);
    
    // Allocate temporary storage
    cudaMalloc(&d_tempStorage, tempStorageBytes);
    
    // Run sum reduction
    cub::DeviceReduce::Sum(d_tempStorage, tempStorageBytes, 
                           d_input, d_output, N);
    
    cudaFree(d_tempStorage);
}

8.6 性能分析与建模

计算强度：

归约的计算强度非常低：

计算量：N-1次加法
内存传输：读取N个元素，写入1个结果
计算强度 ≈ 1 FLOP / 4 bytes = 0.25 FLOPs/byte

这远低于GPU达到峰值性能所需的计算强度（通常>10 FLOPs/byte）。因此，归约是典型的内存带宽受限操作。

理论峰值带宽利用率：

理想的归约应该达到接近理论峰值的内存带宽。实际中，由于同步开销、非合并访问等因素，利用率通常在60-80%之间。

同步开销分析：

每次__syncthreads()引入的开销约为几百个时钟周期。在标准归约中，有log₂(blockSize)次同步。对于blockSize=256，有8次同步。展开最后warp可以将同步次数减少到5次。

Occupancy影响：

归约kernel通常具有较高的occupancy，因为它们使用较少的寄存器和共享内存。高occupancy有助于隐藏内存延迟。

8.7 归约变体

ArgMax/ArgMin：

同时返回最大值/最小值及其索引。需要将值和索引打包为结构体，并自定义归约运算符。

cpp

struct ValIdx {
    float value;
    int index;
};

__device__ ValIdx maxReduce(ValIdx a, ValIdx b) {
    return (a.value >= b.value) ? a : b;
}

直方图：

直方图可以看作是一种特殊的归约，其中bin索引决定数据归属。由于多个线程可能更新同一个bin，需要使用原子操作或共享内存暂存。

稀疏归约：

只对满足条件的元素进行归约。可以在加载阶段添加过滤条件。

三、程序员的视角

3.1 在实际项目中的应用

统计分析：求和、平均值、方差、标准差等统计量都依赖归约。在大数据分析和科学数据处理中，这些操作频繁出现。

深度学习：Batch Normalization需要计算batch内的均值和方差；损失函数的聚合需要归约；梯度同步在多GPU训练中是关键步骤。

物理模拟：总能量、总动量、质心位置等守恒量的计算需要全局归约。在粒子模拟中，边界框的确定需要最大值/最小值归约。

数据库加速：SQL的SUM、COUNT、MAX、MIN等聚合函数本质上是归约操作。GPU加速的数据库引擎需要高效的归约实现。

3.2 常见陷阱与误区

忘记初始化输出：使用atomicAdd前必须将输出初始化为0，否则结果错误。

整数溢出：对于大规模数据的求和，32位整数可能溢出。应使用64位整数或双精度浮点数。

浮点精度损失：不同归约顺序导致不同的舍入误差累积。对于精度敏感的应用，考虑使用Kahan求和或双精度。

Bank Conflict：交错寻址在早期步长中可能导致bank conflict。使用相邻寻址或填充可以避免。

过度使用原子操作：对于大规模数据，直接使用原子操作会导致严重的竞争和序列化。应使用多级归约。

3.3 与其他层级的联系

与前缀和的联系：归约是前缀和的特例（只取最后一个元素）。前缀和的优化技术（如Blelloch算法）也可以应用于归约。

与排序的联系：某些排序算法（如基数排序）内部使用归约来计算直方图。

与稀疏矩阵的联系：稀疏矩阵向量乘法（SpMV）的结果聚合可以看作是一种不规则的归约。

四、代码实战

4.1 完整的归约优化对比

cpp

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

#define N 10000000
#define BLOCK_SIZE 256

// Version 1: Basic interleaved addressing
__global__ void reductionV1(const float* input, float* output, int N) {
    extern __shared__ float shared[];
    int tid = threadIdx.x;
    int idx = blockIdx.x * blockDim.x + threadIdx.x;
    
    shared[tid] = (idx < N) ? input[idx] : 0.0f;
    __syncthreads();
    
    for (int stride = blockDim.x / 2; stride > 0; stride >>= 1) {
        if (tid < stride) {
            shared[tid] += shared[tid + stride];
        }
        __syncthreads();
    }
    
    if (tid == 0) atomicAdd(output, shared[0]);
}

// Version 2: Sequential addressing
__global__ void reductionV2(const float* input, float* output, int N) {
    extern __shared__ float shared[];
    int tid = threadIdx.x;
    int idx = blockIdx.x * (blockDim.x * 2) + threadIdx.x;
    
    float val = (idx < N) ? input[idx] : 0.0f;
    if (idx + blockDim.x < N) val += input[idx + blockDim.x];
    shared[tid] = val;
    __syncthreads();
    
    for (int stride = blockDim.x / 2; stride > 0; stride >>= 1) {
        if (tid < stride) {
            shared[tid] += shared[tid + stride];
        }
        __syncthreads();
    }
    
    if (tid == 0) atomicAdd(output, shared[0]);
}

// Version 3: Unrolled last warp
__global__ void reductionV3(const float* input, float* output, int N) {
    extern __shared__ float shared[];
    int tid = threadIdx.x;
    int idx = blockIdx.x * (blockDim.x * 2) + threadIdx.x;
    
    float val = (idx < N) ? input[idx] : 0.0f;
    if (idx + blockDim.x < N) val += input[idx + blockDim.x];
    shared[tid] = val;
    __syncthreads();
    
    for (int stride = blockDim.x / 2; stride > 32; stride >>= 1) {
        if (tid < stride) {
            shared[tid] += shared[tid + stride];
        }
        __syncthreads();
    }
    
    if (tid < 32) {
        volatile float* vs = shared;
        vs[tid] += vs[tid + 32];
        vs[tid] += vs[tid + 16];
        vs[tid] += vs[tid + 8];
        vs[tid] += vs[tid + 4];
        vs[tid] += vs[tid + 2];
        vs[tid] += vs[tid + 1];
    }
    
    if (tid == 0) atomicAdd(output, shared[0]);
}

// Reference CPU sum
float cpuSum(const float* data, int N) {
    double sum = 0.0;
    for (int i = 0; i < N; i++) sum += data[i];
    return (float)sum;
}

void benchmarkReductions() {
    printf("=== Reduction Optimization Benchmark ===\n");
    
    size_t size = N * sizeof(float);
    float *h_data = (float*)malloc(size);
    for (int i = 0; i < N; i++) h_data[i] = 1.0f;
    
    float refSum = cpuSum(h_data, N);
    printf("Reference sum: %.0f\n", refSum);
    
    float *d_data, *d_result;
    cudaMalloc((void**)&d_data, size);
    cudaMalloc((void**)&d_result, sizeof(float));
    cudaMemcpy(d_data, h_data, size, cudaMemcpyHostToDevice);
    
    int blockSize = BLOCK_SIZE;
    int gridSize = (N + blockSize * 2 - 1) / (blockSize * 2);
    
    cudaEvent_t start, stop;
    cudaEventCreate(&start);
    cudaEventCreate(&stop);
    
    // V1
    cudaMemset(d_result, 0, sizeof(float));
    cudaEventRecord(start);
    for (int i = 0; i < 30; i++) {
        reductionV1<<<gridSize, blockSize, blockSize * sizeof(float)>>>(
            d_data, d_result, N);
    }
    cudaEventRecord(stop);
    cudaEventSynchronize(stop);
    float ms1;
    cudaEventElapsedTime(&ms1, start, stop);
    
    // V2
    cudaMemset(d_result, 0, sizeof(float));
    cudaEventRecord(start);
    for (int i = 0; i < 30; i++) {
        reductionV2<<<gridSize, blockSize, blockSize * sizeof(float)>>>(
            d_data, d_result, N);
    }
    cudaEventRecord(stop);
    cudaEventSynchronize(stop);
    float ms2;
    cudaEventElapsedTime(&ms2, start, stop);
    
    // V3
    cudaMemset(d_result, 0, sizeof(float));
    cudaEventRecord(start);
    for (int i = 0; i < 30; i++) {
        reductionV3<<<gridSize, blockSize, blockSize * sizeof(float)>>>(
            d_data, d_result, N);
    }
    cudaEventRecord(stop);
    cudaEventSynchronize(stop);
    float ms3;
    cudaEventElapsedTime(&ms3, start, stop);
    
    printf("V1 (Interleaved): %.3f ms\n", ms1 / 30);
    printf("V2 (Sequential):  %.3f ms\n", ms2 / 30);
    printf("V3 (Unrolled):    %.3f ms\n", ms3 / 30);
    printf("V2/V1 Speedup:    %.2fx\n", ms1 / ms2);
    printf("V3/V1 Speedup:    %.2fx\n", ms1 / ms3);
    
    cudaFree(d_data);
    cudaFree(d_result);
    free(h_data);
    cudaEventDestroy(start);
    cudaEventDestroy(stop);
}

int main() {
    benchmarkReductions();
    return 0;
}

4.2 ArgMax实现

cpp

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

#define N 1000000
#define BLOCK_SIZE 256

struct ValIdx {
    float value;
    int index;
};

__device__ ValIdx maxReduce(ValIdx a, ValIdx b) {
    return (a.value >= b.value) ? a : b;
}

__global__ void argMaxKernel(const float* input, ValIdx* output, int N) {
    extern __shared__ ValIdx shared[];
    
    int tid = threadIdx.x;
    int idx = blockIdx.x * blockDim.x + threadIdx.x;
    
    ValIdx local;
    if (idx < N) {
        local.value = input[idx];
        local.index = idx;
    } else {
        local.value = -FLT_MAX;
        local.index = -1;
    }
    shared[tid] = local;
    __syncthreads();
    
    for (int stride = blockDim.x / 2; stride > 32; stride >>= 1) {
        if (tid < stride) {
            shared[tid] = maxReduce(shared[tid], shared[tid + stride]);
        }
        __syncthreads();
    }
    
    if (tid < 32) {
        volatile ValIdx* vs = shared;
        vs[tid] = maxReduce(vs[tid], vs[tid + 32]);
        vs[tid] = maxReduce(vs[tid], vs[tid + 16]);
        vs[tid] = maxReduce(vs[tid], vs[tid + 8]);
        vs[tid] = maxReduce(vs[tid], vs[tid + 4]);
        vs[tid] = maxReduce(vs[tid], vs[tid + 2]);
        vs[tid] = maxReduce(vs[tid], vs[tid + 1]);
    }
    
    if (tid == 0) {
        // Use atomicCAS to find global max
        ValIdx* globalMax = output;
        ValIdx old = *globalMax;
        ValIdx assumed;
        do {
            assumed = old;
            if (shared[0].value <= assumed.value) break;
            old = atomicCAS((unsigned long long*)globalMax, 
                           *(unsigned long long*)&assumed,
                           *(unsigned long long*)&shared[0]);
        } while (assumed.value == old.value && assumed.index == old.index);
    }
}

int main() {
    float *h_data = (float*)malloc(N * sizeof(float));
    int maxIdx = N / 3;
    for (int i = 0; i < N; i++) {
        h_data[i] = (float)i;
    }
    
    float *d_data;
    ValIdx *d_result;
    cudaMalloc((void**)&d_data, N * sizeof(float));
    cudaMalloc((void**)&d_result, sizeof(ValIdx));
    
    cudaMemcpy(d_data, h_data, N * sizeof(float), cudaMemcpyHostToDevice);
    
    // Initialize result to -inf
    ValIdx h_init = {-FLT_MAX, -1};
    cudaMemcpy(d_result, &h_init, sizeof(ValIdx), cudaMemcpyHostToDevice);
    
    int blockSize = BLOCK_SIZE;
    int gridSize = (N + blockSize - 1) / blockSize;
    
    argMaxKernel<<<gridSize, blockSize, blockSize * sizeof(ValIdx)>>>(
        d_data, d_result, N);
    
    ValIdx h_result;
    cudaMemcpy(&h_result, d_result, sizeof(ValIdx), cudaMemcpyDeviceToHost);
    
    printf("ArgMax: index=%d, value=%.0f\n", h_result.index, h_result.value);
    printf("Expected: index=%d, value=%.0f\n", N-1, (float)(N-1));
    
    cudaFree(d_data);
    cudaFree(d_result);
    free(h_data);
    
    return 0;
}

五、深入思考

5.1 历史背景

并行归约的研究可以追溯到早期的并行计算机时代。Hillis和Steele在1986年提出了前缀和的并行算法，为归约奠定了基础。在GPU领域，Mark Harris的NVIDIA技术报告"Parallel Prefix Sum (Scan) with CUDA"是经典文献，系统地介绍了归约和扫描的优化技术。

随着GPU架构的演进，归约的实现也在不断优化。Fermi引入的原子操作简化了多Block归约；Volta的独立线程调度使得warp级同步更加灵活；Ampere加速的FP32原子操作改变了归约的设计权衡。

5.2 前沿方向

硬件加速归约：一些新兴的GPU架构正在探索硬件级别的归约支持，如专用的归约指令或网络。这将进一步降低归约的延迟。

分布式归约：在多GPU和多节点系统中，归约需要跨设备协调。NVLink和InfiniBand的高速互联使得分布式归约更加高效。

近似归约：对于某些应用（如机器学习），近似的归约结果已经足够。随机采样和Sketching技术可以在保证精度的同时大幅减少通信量。

5.3 延伸阅读

Mark Harris, "Optimizing Parallel Reduction in CUDA", NVIDIA Developer Technology

CUB Library Documentation - DeviceReduce

Blelloch, "Prefix Sums and Their Applications", 1990

NVIDIA CUDA Samples - reduction

六、实验笔记

实验1：归约优化效果量化

目的：测量每种优化技术的贡献。

步骤：

实现V1-V3三个版本的归约

使用Nsight Compute分析每个版本的同步次数、内存事务和占用率

绘制加速比对比图

预期：V2相比V1应有1.2-1.5倍加速，主要来自更好的内存访问模式。V3相比V2应有进一步的1.1-1.3倍加速，来自减少的同步开销。

实验2：原子操作vs树形归约

目的：比较原子操作和树形归约在不同数据规模下的性能。

步骤：

实现纯原子操作版本和多级归约版本

测试N = 10^3, 10^4, ..., 10^8的性能

找出原子操作更优的临界点

预期：对于小数据（N < 10^4），原子操作可能更快，因为开销小。对于大数据，树形归约明显优于原子操作。

七、本章小结

归约是并行计算中最基础的模式之一，它将大量数据通过二元运算符合并为单个值。本章从归约的数学定义出发，讨论了结合律和交换律对并行化的重要性。

树形归约是并行归约的核心算法，通过将数据组织为二叉树结构，将并行深度从O(N)降为O(log N)。基于共享内存的实现利用了GPU的快速片上存储器，显著减少了全局内存访问。

我们详细介绍了多种优化技术：相邻寻址改善了内存访问模式，展开最后warp消除了冗余同步，多元素处理减少了kernel启动开销，原子操作提供了简化的替代方案。通过这些优化，归约性能可以接近硬件的理论峰值带宽。

多Block归约需要协调多个Block的部分结果，常用的策略包括两阶段归约和使用CUB库。在生产代码中，推荐使用经过充分测试的库函数而非手写实现。

归约的变体如ArgMax、直方图等在实际应用中同样重要。通过自定义归约运算符，可以将归约模式扩展到各种场景。

八、习题

结合律验证：浮点加法是否严格满足结合律？给出一个反例说明(a + b) + c ≠ a + (b + c)的情况。

归约树深度：对于blockSize=512，归约树有多少层？如果展开最后warp，还剩多少次同步？

Bank Conflict分析：在交错寻址的归约中，第一次迭代（stride=blockSize/2）是否存在bank conflict？为什么？

原子操作竞争：如果有10000个线程同时执行atomicAdd(&globalSum, value)，估计序列化导致的性能下降。

ArgMin实现：修改ArgMax kernel实现ArgMin，同时返回最小值和其索引。

多元素负载均衡：在多元素处理的归约中，如果N不能被gridSize整除，如何确保负载均衡？

双精度归约：将归约kernel改为支持double类型。需要注意哪些变化？

二维归约：对于一个M×N的矩阵，如何高效地计算所有元素的和？考虑行优先和列优先的访问模式。

条件归约：实现一个归约，只对大于阈值的元素求和。如何在加载阶段过滤数据？

综合优化：设计一个针对10^9元素的FP32数组的归约方案。考虑Grid大小、Block大小、多阶段策略等因素，最大化吞吐量。

第8章 并行模式：归约 - 从数据到结果 (Parallel Patterns: Reduction - From Data to Result)