第5章线程同步 - 协调并行线程（Thread Synchronization）

一、导读

1.1 本章学习目标

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

理解线程同步在并行计算中的必要性和挑战

掌握 __syncthreads() 的工作原理和使用方法

理解屏障同步（Barrier Synchronization）的概念和实现机制

学会使用共享内存配合线程同步实现高效的数据交换

掌握原子操作（Atomic Operations）的类型和应用场景

理解竞争条件（Race Condition）的成因和解决方法

能够识别和避免死锁（Deadlock）问题

掌握 warp-level 同步原语的使用方法

理解 CUDA 内存模型与同步的关系

能够设计高效的并行归约和扫描算法

1.2 为什么本章如此重要

在并行计算中，线程同步是确保正确性的关键机制。当多个线程同时访问共享资源时，如果没有适当的同步，可能会导致数据竞争、结果不一致等严重问题。

CUDA 提供了多层次的同步机制：

块内同步：__syncthreads() 用于同一线程块内的线程同步
原子操作：用于全局范围内的安全内存更新
Warp-level 原语：用于 warp 内的高效协作

理解这些同步机制的工作原理和适用场景，是编写正确且高效的 CUDA 程序的基础。错误的同步可能导致程序崩溃、结果错误或性能下降；而过度的同步则可能限制并行度，降低效率。

本章还将介绍一些经典的并行算法模式，如归约（Reduction）、扫描（Scan）和前缀和（Prefix Sum），这些模式广泛应用于各种并行计算任务中。

1.3 核心问题引导

在深入学习之前，请思考以下问题：

为什么线程需要同步？不同步会发生什么？

__syncthreads() 的作用范围是什么？它能跨线程块同步吗？

什么是竞争条件？如何检测和避免？

原子操作的性能代价有多大？什么时候应该使用它们？

什么是死锁？在 CUDA 中如何避免死锁？

warp-level 同步与块级同步有什么区别？

如何使用同步机制实现高效的并行归约？

共享内存和同步如何配合使用以实现数据复用？

二、核心概念详解

5.1 线程同步的必要性

并行执行的不确定性

在 GPU 上，线程的执行顺序是不确定的。即使是同一个内核，每次运行时线程的调度顺序也可能不同。这种不确定性带来了两个主要问题：

数据依赖问题：如果一个线程需要读取另一个线程写入的数据，必须确保写操作在读操作之前完成。

竞争条件：当多个线程同时读写同一内存位置时，最终结果取决于执行顺序，这可能导致不可预测的行为。

示例：没有同步的问题

cpp

__global__ void unsafeKernel(float *sharedData, float *input) {
    __shared__ float temp[256];
    int tid = threadIdx.x;
    
    // 线程 A 写入
    temp[tid] = input[tid];
    
    // 线程 B 读取 - 可能读到未初始化的数据！
    if (tid > 0) {
        temp[tid] += temp[tid - 1]; // 危险！temp[tid-1] 可能还未写入
    }
}

在这个例子中，线程 tid 依赖于线程 tid-1 的结果，但由于没有同步，temp[tid-1] 可能尚未被写入。

5.2 __syncthreads() 屏障同步

基本用法

__syncthreads() 是 CUDA 中最基本的同步原语。它在当前线程块内建立一个屏障，所有线程必须到达这个屏障才能继续执行。

cpp

__global__ void synchronizedKernel(float *output, float *input, int n) {
    __shared__ float sharedData[256];
    int tid = threadIdx.x;
    int globalIdx = blockIdx.x * blockDim.x + threadIdx.x;
    
    // 阶段 1：加载数据到共享内存
    if (globalIdx < n) {
        sharedData[tid] = input[globalIdx];
    }
    
    // 屏障：等待所有线程完成加载
    __syncthreads();
    
    // 阶段 2：处理共享内存中的数据
    if (tid < blockDim.x / 2) {
        sharedData[tid] += sharedData[tid + blockDim.x / 2];
    }
    
    // 再次同步
    __syncthreads();
    
    // 阶段 3：写回结果
    if (tid == 0) {
        output[blockIdx.x] = sharedData[0];
    }
}

工作原理

当线程执行到 __syncthreads() 时：

该线程暂停执行

硬件跟踪每个线程块中已到达屏障的线程数

当块内所有活动线程都到达屏障时，所有线程同时释放

线程继续执行屏障后的代码

重要特性

块内作用：只同步同一线程块内的线程，不同块之间不同步
隐式内存屏障：__syncthreads() 也充当内存屏障，确保之前的内存操作对后续操作可见
条件同步的危险：如果部分线程跳过 __syncthreads()，会导致死锁

条件同步的陷阱

cpp

// 错误示例：可能导致死锁
__global__ void badConditionalSync(float *data) {
    __shared__ float temp[256];
    
    if (threadIdx.x < 128) {
        temp[threadIdx.x] = data[threadIdx.x];
        __syncthreads(); // 只有前 128 个线程执行到这里
        // 后 128 个线程永远不会到达这里 → 死锁！
    } else {
        // 这些线程跳过了上面的 __syncthreads()
        __syncthreads(); // 等待永远无法完成的同步
    }
}

// 正确示例：所有线程都执行同步
__global__ void goodConditionalSync(float *data) {
    __shared__ float temp[256];
    
    if (threadIdx.x < 128) {
        temp[threadIdx.x] = data[threadIdx.x];
    }
    __syncthreads(); // 所有线程都执行
    
    if (threadIdx.x < 128) {
        // 处理数据
        temp[threadIdx.x] *= 2.0f;
    }
    __syncthreads(); // 所有线程都执行
}

5.3 原子操作（Atomic Operations）

什么是原子操作？

原子操作是不可分割的操作，在执行过程中不会被其他线程中断。CUDA 提供了一系列原子函数，用于在全局或共享内存中安全地更新值。

常用原子函数

函数	描述
`atomicAdd(addr, val)`	原子加法：*addr += val
`atomicSub(addr, val)`	原子减法
`atomicExch(addr, val)`	原子交换
`atomicMin(addr, val)`	原子最小值
`atomicMax(addr, val)`	原子最大值
`atomicAnd(addr, val)`	原子按位与
`atomicOr(addr, val)`	原子按位或
`atomicXor(addr, val)`	原子按位异或
`atomicCAS(addr, compare, val)`	原子比较并交换

使用示例

cpp

__global__ void atomicHistogram(int *histogram, int *data, int n, int numBins) {
    int idx = blockIdx.x * blockDim.x + threadIdx.x;
    if (idx < n) {
        int bin = data[idx] % numBins;
        atomicAdd(&histogram[bin], 1); // 安全地增加计数
    }
}

性能考虑

原子操作的性能代价：

全局内存原子操作：较慢，因为需要独占访问全局内存
共享内存原子操作：较快，但仍比非原子操作慢
冲突程度：越多线程同时访问同一地址，性能下降越明显

优化策略

使用共享内存暂存：先在共享内存中累积，最后再原子更新全局内存

减少冲突：使用私有计数器或分散访问模式

选择合适的数据类型：某些架构对 32 位原子操作优化更好

cpp

// 优化的直方图计算
__global__ void optimizedHistogram(int *histogram, int *data, int n, int numBins) {
    extern __shared__ int sharedHist[];
    
    int tid = threadIdx.x;
    int idx = blockIdx.x * blockDim.x + threadIdx.x;
    
    // 初始化共享内存直方图
    for (int i = tid; i < numBins; i += blockDim.x) {
        sharedHist[i] = 0;
    }
    __syncthreads();
    
    // 在共享内存中累积
    if (idx < n) {
        int bin = data[idx] % numBins;
        atomicAdd(&sharedHist[bin], 1);
    }
    __syncthreads();
    
    // 合并到全局直方图
    if (tid < numBins) {
        atomicAdd(&histogram[tid], sharedHist[tid]);
    }
}

5.4 竞争条件（Race Condition）

定义

竞争条件发生在多个线程并发访问同一数据，且至少有一个是写操作时。最终结果取决于线程的执行顺序。

常见场景

累加操作：多个线程同时更新同一计数器

链表操作：多个线程同时修改链表指针

标志位检查：检查和设置标志位不是原子操作

检测方法

CUDA-Memcheck：NVIDIA 提供的工具可以检测数据竞争
重复运行：竞争条件可能导致间歇性错误
代码审查：仔细分析共享数据的访问模式

解决方法

使用原子操作：适用于简单的读-改-写操作

使用同步：确保操作的顺序性

避免共享：为每个线程分配独立的数据空间

使用归约模式：将全局操作分解为局部操作

5.5 死锁（Deadlock）

什么是死锁？

死锁发生在两个或多个线程互相等待对方释放资源时，导致所有线程都无法继续执行。

CUDA 中的死锁原因

条件同步：部分线程跳过 __syncthreads()

嵌套同步：在共享内存不足时，过多的同步点可能导致资源耗尽

循环依赖：线程 A 等待线程 B，线程 B 等待线程 A

避免死锁的最佳实践

确保所有线程执行相同的同步路径

避免在条件语句中使用同步

保持同步逻辑简单明了

测试边界条件

5.6 Warp-Level 同步原语

Volta 架构引入了独立的线程调度，同时也提供了 warp-level 的同步原语。

__syncwarp()

cpp

__global__ void warpLevelSync(float *data) {
    int tid = threadIdx.x;
    int warpId = tid / 32;
    int laneId = tid % 32;
    
    // warp 内同步
    data[tid] = ...;
    __syncwarp(); // 只同步当前 warp 内的线程
    
    // 使用 warp shuffle 进行数据交换
    float neighbor = __shfl_down_sync(0xFFFFFFFF, data[tid], 1);
}

Warp Shuffle 函数

Warp shuffle 允许 warp 内的线程直接交换寄存器数据，无需通过共享内存。

函数	描述
`__shfl_sync(mask, var, srcLane)`	从指定线程获取值
`__shfl_up_sync(mask, var, delta)`	从索引较小的线程获取值
`__shfl_down_sync(mask, var, delta)`	从索引较大的线程获取值
`__shfl_xor_sync(mask, var, laneMask)`	从 XOR 掩码计算的线程获取值

使用 warp shuffle 实现归约

cpp

__device__ float warpReduce(float val) {
    for (int offset = 16; offset > 0; offset /= 2) {
        val += __shfl_down_sync(0xFFFFFFFF, val, offset);
    }
    return val;
}

__global__ void warpBasedReduce(float *input, float *output, int n) {
    __shared__ float sharedData[32]; // 每个 warp 一个槽位
    
    int tid = threadIdx.x;
    int globalIdx = blockIdx.x * blockDim.x + threadIdx.x;
    int warpId = tid / 32;
    int laneId = tid % 32;
    
    float val = (globalIdx < n) ? input[globalIdx] : 0.0f;
    
    // Warp 内归约
    val = warpReduce(val);
    
    // 每个 warp 的第一个线程写入共享内存
    if (laneId == 0) {
        sharedData[warpId] = val;
    }
    __syncthreads();
    
    // 最后一个 warp 进行最终归约
    if (warpId == 0) {
        val = (tid < blockDim.x / 32) ? sharedData[laneId] : 0.0f;
        val = warpReduce(val);
        if (tid == 0) {
            output[blockIdx.x] = val;
        }
    }
}

5.7 CUDA 内存模型与同步

内存可见性

在 CUDA 中，一个线程写入的数据何时对其他线程可见？这取决于内存类型和同步机制：

寄存器：线程私有，永远不可见

共享内存：同一块内的线程可见，需要 __syncthreads() 确保可见性

全局内存：所有线程可见，但可能需要内存屏障

常量/纹理内存：只读，无需考虑可见性

内存屏障

__syncthreads() 不仅同步线程执行，还充当内存屏障：

确保屏障前的所有内存操作在屏障后可见
防止编译器和硬件的重排序优化破坏程序语义

volatile 关键字

对于频繁更新的共享变量，可以使用 volatile 告诉编译器不要缓存该变量：

cpp

__shared__ volatile int flag;

三、程序员的视角

5.1 在实际项目中的应用

并行归约

归约是将多个值合并为一个值的操作，如求和、求最大值等。这是最常见的并行模式之一。

并行扫描（Prefix Sum）

扫描计算数组的前缀和，广泛应用于排序、压缩等算法。

粒子系统更新

在物理模拟中，粒子之间的相互作用需要同步机制来确保一致性。

5.2 常见陷阱与误区

误区一：__syncthreads() 可以跨块同步

事实：__syncthreads() 只在块内有效。如果需要跨块同步，必须分多个内核调用。

误区二：原子操作总是安全的

事实：原子操作只能保证单个操作的原子性，不能保证多个操作的组合原子性。

陷阱：在条件分支中使用同步

如果控制流导致部分线程跳过同步点，会导致死锁。

陷阱：过度使用原子操作

频繁的原子操作会序列化执行，抵消并行带来的收益。

5.3 与其他层级的联系

第3章：线程模型决定了同步的范围和方式
第4章：共享内存的使用通常配合同步机制
第8章：同步是影响性能的重要因素

四、代码实战

5.1 完整的并行归约实现

cpp

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

// 基础归约内核
__global__ void reduceBasic(float *input, float *output, int n) {
    extern __shared__ float sharedData[];
    
    int tid = threadIdx.x;
    int globalIdx = blockIdx.x * blockDim.x + threadIdx.x;
    
    // 加载数据
    sharedData[tid] = (globalIdx < n) ? input[globalIdx] : 0.0f;
    __syncthreads();
    
    // 树形归约
    for (int stride = blockDim.x / 2; stride > 0; stride >>= 1) {
        if (tid < stride) {
            sharedData[tid] += sharedData[tid + stride];
        }
        __syncthreads();
    }
    
    // 写入结果
    if (tid == 0) {
        output[blockIdx.x] = sharedData[0];
    }
}

// 优化的归约：减少同步次数
__global__ void reduceOptimized(float *input, float *output, int n) {
    extern __shared__ float sharedData[];
    
    int tid = threadIdx.x;
    int globalIdx = blockIdx.x * blockDim.x + threadIdx.x;
    
    // 加载数据并在寄存器中进行初步归约
    float val = (globalIdx < n) ? input[globalIdx] : 0.0f;
    if (globalIdx + blockDim.x < n) {
        val += input[globalIdx + blockDim.x];
    }
    
    sharedData[tid] = val;
    __syncthreads();
    
    // 树形归约
    for (int stride = blockDim.x / 2; stride > 0; stride >>= 1) {
        if (tid < stride) {
            sharedData[tid] += sharedData[tid + stride];
        }
        __syncthreads();
    }
    
    if (tid == 0) {
        output[blockIdx.x] = sharedData[0];
    }
}

// 使用 warp shuffle 的归约
__device__ float warpReduceSum(float val) {
    for (int offset = 16; offset > 0; offset /= 2) {
        val += __shfl_down_sync(0xFFFFFFFF, val, offset);
    }
    return val;
}

__global__ void reduceWarpShuffle(float *input, float *output, int n) {
    __shared__ float sharedData[32];
    
    int tid = threadIdx.x;
    int globalIdx = blockIdx.x * blockDim.x + threadIdx.x;
    int warpId = tid / 32;
    int laneId = tid % 32;
    
    float val = (globalIdx < n) ? input[globalIdx] : 0.0f;
    val = warpReduceSum(val);
    
    if (laneId == 0) {
        sharedData[warpId] = val;
    }
    __syncthreads();
    
    if (warpId == 0) {
        val = (tid < blockDim.x / 32) ? sharedData[laneId] : 0.0f;
        val = warpReduceSum(val);
        if (tid == 0) {
            output[blockIdx.x] = val;
        }
    }
}

#define CUDA_CHECK(call) \
    do { \
        cudaError_t err = call; \
        if (err != cudaSuccess) { \
            printf("CUDA error: %s\n", cudaGetErrorString(err)); \
            exit(1); \
        } \
    } while(0)

int main() {
    int n = 1 << 20;
    size_t size = n * sizeof(float);
    
    float *h_input = (float *)malloc(size);
    float expectedSum = 0.0f;
    
    for (int i = 0; i < n; i++) {
        h_input[i] = 1.0f;
        expectedSum += 1.0f;
    }
    
    float *d_input;
    CUDA_CHECK(cudaMalloc((void **)&d_input, size));
    CUDA_CHECK(cudaMemcpy(d_input, h_input, size, cudaMemcpyHostToDevice));
    
    int threadsPerBlock = 256;
    int blocksPerGrid = (n + threadsPerBlock - 1) / threadsPerBlock;
    size_t outputSize = blocksPerGrid * sizeof(float);
    
    float *d_output;
    CUDA_CHECK(cudaMalloc((void **)&d_output, outputSize));
    
    cudaEvent_t start, stop;
    CUDA_CHECK(cudaEventCreate(&start));
    CUDA_CHECK(cudaEventCreate(&stop));
    
    // 测试基础归约
    int sharedMemSize = threadsPerBlock * sizeof(float);
    CUDA_CHECK(cudaEventRecord(start));
    reduceBasic<<<blocksPerGrid, threadsPerBlock, sharedMemSize>>>(d_input, d_output, n);
    CUDA_CHECK(cudaEventRecord(stop));
    CUDA_CHECK(cudaEventSynchronize(stop));
    
    float basicTime = 0;
    CUDA_CHECK(cudaEventElapsedTime(&basicTime, start, stop));
    printf("Basic reduce time: %.3f ms\n", basicTime);
    
    // 测试优化归约
    CUDA_CHECK(cudaEventRecord(start));
    reduceOptimized<<<blocksPerGrid, threadsPerBlock, sharedMemSize>>>(d_input, d_output, n);
    CUDA_CHECK(cudaEventRecord(stop));
    CUDA_CHECK(cudaEventSynchronize(stop));
    
    float optimizedTime = 0;
    CUDA_CHECK(cudaEventElapsedTime(&optimizedTime, start, stop));
    printf("Optimized reduce time: %.3f ms\n", optimizedTime);
    
    // 测试 warp shuffle 归约
    CUDA_CHECK(cudaEventRecord(start));
    reduceWarpShuffle<<<blocksPerGrid, threadsPerBlock, sharedMemSize>>>(d_input, d_output, n);
    CUDA_CHECK(cudaEventRecord(stop));
    CUDA_CHECK(cudaEventSynchronize(stop));
    
    float warpTime = 0;
    CUDA_CHECK(cudaEventElapsedTime(&warpTime, start, stop));
    printf("Warp shuffle reduce time: %.3f ms\n", warpTime);
    
    // 验证结果
    float *h_output = (float *)malloc(outputSize);
    CUDA_CHECK(cudaMemcpy(h_output, d_output, outputSize, cudaMemcpyDeviceToHost));
    
    float sum = 0.0f;
    for (int i = 0; i < blocksPerGrid; i++) {
        sum += h_output[i];
    }
    
    printf("Sum: %.0f (expected: %.0f)\n", sum, expectedSum);
    printf("Verification: %s\n", (fabs(sum - expectedSum) < 1e-3) ? "PASSED" : "FAILED");
    
    CUDA_CHECK(cudaFree(d_input));
    CUDA_CHECK(cudaFree(d_output));
    CUDA_CHECK(cudaEventDestroy(start));
    CUDA_CHECK(cudaEventDestroy(stop));
    free(h_input);
    free(h_output);
    
    return 0;
}

5.2 原子操作直方图

cpp

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

// 简单原子直方图
__global__ void simpleHistogram(int *histogram, int *data, int n, int numBins) {
    int idx = blockIdx.x * blockDim.x + threadIdx.x;
    if (idx < n) {
        int bin = data[idx] % numBins;
        atomicAdd(&histogram[bin], 1);
    }
}

// 优化的直方图：使用共享内存
__global__ void optimizedHistogram(int *histogram, int *data, int n, int numBins) {
    extern __shared__ int sharedHist[];
    
    int tid = threadIdx.x;
    int idx = blockIdx.x * blockDim.x + threadIdx.x;
    
    // 初始化共享内存
    for (int i = tid; i < numBins; i += blockDim.x) {
        sharedHist[i] = 0;
    }
    __syncthreads();
    
    // 累积到共享内存
    if (idx < n) {
        int bin = data[idx] % numBins;
        atomicAdd(&sharedHist[bin], 1);
    }
    __syncthreads();
    
    // 合并到全局内存
    if (tid < numBins) {
        atomicAdd(&histogram[tid], sharedHist[tid]);
    }
}

#define CUDA_CHECK(call) \
    do { \
        cudaError_t err = call; \
        if (err != cudaSuccess) { \
            printf("CUDA error: %s\n", cudaGetErrorString(err)); \
            exit(1); \
        } \
    } while(0)

int main() {
    int n = 1 << 20;
    int numBins = 256;
    size_t dataSize = n * sizeof(int);
    size_t histSize = numBins * sizeof(int);
    
    int *h_data = (int *)malloc(dataSize);
    int *h_histogram = (int *)calloc(numBins, sizeof(int));
    int *h_histogramOpt = (int *)calloc(numBins, sizeof(int));
    
    for (int i = 0; i < n; i++) {
        h_data[i] = rand() % 1000;
    }
    
    // CPU 参考结果
    for (int i = 0; i < n; i++) {
        h_histogram[h_data[i] % numBins]++;
    }
    
    int *d_data, *d_histogram, *d_histogramOpt;
    CUDA_CHECK(cudaMalloc((void **)&d_data, dataSize));
    CUDA_CHECK(cudaMalloc((void **)&d_histogram, histSize));
    CUDA_CHECK(cudaMalloc((void **)&d_histogramOpt, histSize));
    
    CUDA_CHECK(cudaMemcpy(d_data, h_data, dataSize, cudaMemcpyHostToDevice));
    
    int threadsPerBlock = 256;
    int blocksPerGrid = (n + threadsPerBlock - 1) / threadsPerBlock;
    
    // 测试简单版本
    CUDA_CHECK(cudaMemset(d_histogram, 0, histSize));
    simpleHistogram<<<blocksPerGrid, threadsPerBlock>>>(d_histogram, d_data, n, numBins);
    
    float *h_result = (float *)malloc(histSize);
    CUDA_CHECK(cudaMemcpy(h_result, d_histogram, histSize, cudaMemcpyDeviceToHost));
    
    bool correct = true;
    for (int i = 0; i < numBins; i++) {
        if (h_result[i] != h_histogram[i]) {
            correct = false;
            break;
        }
    }
    printf("Simple histogram: %s\n", correct ? "PASSED" : "FAILED");
    
    // 测试优化版本
    int sharedMemSize = numBins * sizeof(int);
    CUDA_CHECK(cudaMemset(d_histogramOpt, 0, histSize));
    optimizedHistogram<<<blocksPerGrid, threadsPerBlock, sharedMemSize>>>(
        d_histogramOpt, d_data, n, numBins);
    
    CUDA_CHECK(cudaMemcpy(h_result, d_histogramOpt, histSize, cudaMemcpyDeviceToHost));
    
    correct = true;
    for (int i = 0; i < numBins; i++) {
        if (h_result[i] != h_histogram[i]) {
            correct = false;
            break;
        }
    }
    printf("Optimized histogram: %s\n", correct ? "PASSED" : "FAILED");
    
    CUDA_CHECK(cudaFree(d_data));
    CUDA_CHECK(cudaFree(d_histogram));
    CUDA_CHECK(cudaFree(d_histogramOpt));
    free(h_data);
    free(h_histogram);
    free(h_histogramOpt);
    free(h_result);
    
    return 0;
}

5.3 并行前缀和（Scan）

cpp

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

// Hillis-Steele 扫描算法
__global__ void scanKernel(float *output, const float *input, int n) {
    extern __shared__ float sharedData[];
    
    int tid = threadIdx.x;
    int globalIdx = blockIdx.x * blockDim.x + threadIdx.x;
    
    // 加载数据
    if (globalIdx < n) {
        sharedData[tid] = input[globalIdx];
    } else {
        sharedData[tid] = 0.0f;
    }
    __syncthreads();
    
    // Upsweep 阶段
    for (int stride = 1; stride < blockDim.x; stride *= 2) {
        float temp = 0.0f;
        if (tid >= stride) {
            temp = sharedData[tid - stride];
        }
        __syncthreads();
        if (tid >= stride) {
            sharedData[tid] += temp;
        }
        __syncthreads();
    }
    
    // Downsweep 阶段
    for (int stride = blockDim.x / 2; stride > 0; stride /= 2) {
        float temp = sharedData[tid];
        if (tid >= stride) {
            temp -= sharedData[tid - stride];
        }
        __syncthreads();
        if (tid >= stride) {
            sharedData[tid] = temp;
        }
        __syncthreads();
    }
    
    // 写回结果
    if (globalIdx < n) {
        output[globalIdx] = sharedData[tid];
    }
}

#define CUDA_CHECK(call) \
    do { \
        cudaError_t err = call; \
        if (err != cudaSuccess) { \
            printf("CUDA error: %s\n", cudaGetErrorString(err)); \
            exit(1); \
        } \
    } while(0)

int main() {
    int n = 1024;
    size_t size = n * sizeof(float);
    
    float *h_input = (float *)malloc(size);
    float *h_output = (float *)malloc(size);
    float *h_expected = (float *)malloc(size);
    
    for (int i = 0; i < n; i++) {
        h_input[i] = 1.0f;
    }
    
    // CPU 参考结果
    h_expected[0] = h_input[0];
    for (int i = 1; i < n; i++) {
        h_expected[i] = h_expected[i-1] + h_input[i];
    }
    
    float *d_input, *d_output;
    CUDA_CHECK(cudaMalloc((void **)&d_input, size));
    CUDA_CHECK(cudaMalloc((void **)&d_output, size));
    
    CUDA_CHECK(cudaMemcpy(d_input, h_input, size, cudaMemcpyHostToDevice));
    
    int threadsPerBlock = 1024;
    int blocksPerGrid = 1;
    int sharedMemSize = threadsPerBlock * sizeof(float);
    
    scanKernel<<<blocksPerGrid, threadsPerBlock, sharedMemSize>>>(d_output, d_input, n);
    CUDA_CHECK(cudaDeviceSynchronize());
    
    CUDA_CHECK(cudaMemcpy(h_output, d_output, size, cudaMemcpyDeviceToHost));
    
    // 验证
    bool correct = true;
    for (int i = 0; i < n; i++) {
        if (fabs(h_output[i] - h_expected[i]) > 1e-3) {
            printf("Mismatch at index %d: got %.0f, expected %.0f\n", 
                   i, h_output[i], h_expected[i]);
            correct = false;
            break;
        }
    }
    printf("Scan result: %s\n", correct ? "PASSED" : "FAILED");
    
    CUDA_CHECK(cudaFree(d_input));
    CUDA_CHECK(cudaFree(d_output));
    free(h_input);
    free(h_output);
    free(h_expected);
    
    return 0;
}

五、深入思考

5.1 历史背景

并行同步的研究可以追溯到早期的多处理器系统。Barrier 同步的概念在 1980 年代就被提出，而原子操作则是硬件层面的基本原语。

CUDA 的创新在于将这些概念适配到大规模并行架构中。__syncthreads() 的设计考虑了 GPU 的硬件特性，使其能够在数千个线程之间高效同步。

5.2 前沿方向

当前的研究正在探索：

细粒度同步：更灵活的同步范围，如 warp-level 和 cooperative groups
异步编程模型：减少同步开销的方法
硬件事务内存：自动处理竞争的机制

5.3 延伸阅读

"Parallel Programming in OpenMP" - 同步概念的通用介绍
NVIDIA CUDA C Programming Guide - Synchronization 章节
Cooperative Groups API 文档

六、实验笔记

实验 1：比较不同归约算法的性能

目标：评估基础归约、优化归约和 warp shuffle 归约的性能差异。

步骤：

实现三种归约算法

在不同数据规模下测试

记录执行时间和准确率

分析性能差异的原因

实验 2：探究死锁场景

目标：理解条件同步导致死锁的机制。

步骤：

编写故意造成死锁的内核

观察程序行为

修复死锁并验证

总结避免死锁的规则

七、本章小结

本章详细介绍了 CUDA 中的线程同步机制。我们学习了 __syncthreads() 屏障同步、原子操作、warp-level 原语等多种同步工具。

关键要点包括：

__syncthreads() 用于块内同步，所有线程必须执行相同的同步路径
原子操作用于安全的并发内存更新，但有性能代价
Warp shuffle 提供了高效的 warp 内数据交换机制
竞争条件和死锁是并行编程中的常见问题，需要仔细设计来避免
归约和扫描是重要的并行算法模式

在下一章中，我们将探讨流与并发，学习如何重叠计算和数据传输以进一步提高性能。

八、习题

解释 __syncthreads() 的作用范围和注意事项。

什么是竞争条件？举一个 CUDA 中的例子。

列举常用的原子操作函数及其用途。

为什么在条件分支中使用 __syncthreads() 可能导致死锁？

Warp shuffle 与共享内存相比有什么优势？

编写一个内核，使用原子操作计算数组中大于阈值的元素个数。

解释 Hillis-Steele 扫描算法的工作原理。

如何优化直方图计算以减少原子操作的冲突？

什么是内存屏障？它与线程同步有什么关系？

设计一个实验来测量原子操作与非原子操作的性能差异。

第5章 线程同步 - 协调并行线程（Thread Synchronization）