Lecture 6: 并行模式四 —— 归约：树形归约与 warp 级原语 (对应 Lab 4 / Lab 5: List Reduction)

if (tid < s) w += shared[tid + s];
__syncthreads();

// 文件: reduce_v1_interleaved.cu
// 编译: nvcc -O3 -arch=sm_80 reduce_v1_interleaved.cu -o reduce_v1
// 运行: ./reduce_v1                默认 N = 100000000（1 亿个 float，400 MB）
//       ./reduce_v1 16777216       自定义元素个数 N（任意正整数均可）
//
// 本程序包含两种"交错寻址（interleaved addressing）"归约：
//   变体 1A：Sum25 讲义形式，shared[2*BLOCK_SIZE]，每线程装载 2 个元素，
//            for (stride = 1; stride <= blockDim.x; stride *= 2) if (t % stride == 0)
//   变体 1B：Kirk/Hwu 经典形式，shared[BLOCK_SIZE]，每线程装载 1 个元素，
//            for (s = 1; s < blockDim.x; s *= 2) if (tid % (2*s) == 0)
// 两者逻辑等价、工作量都是 N-1 次加法，但都会造成严重的 warp 发散。

#include <cstdio>
#include <cstdlib>
#include <cmath>
#include <utility>
#include <cuda_runtime.h>

#define BLOCK_SIZE 1024

#define CUDA_CHECK(call)                                                        \
    do {                                                                        \
        cudaError_t err__ = (call);                                             \
        if (err__ != cudaSuccess) {                                             \
            fprintf(stderr, "CUDA error %s:%d: %s\n", __FILE__, __LINE__,       \
                    cudaGetErrorString(err__));                                 \
            exit(EXIT_FAILURE);                                                 \
        }                                                                       \
    } while (0)

// ---------------------------------------------------------------------------
// 变体 1A：讲义 Sum25 第 31 页的交错寻址写法
//   shared 大小 2*BLOCK_SIZE，stride 从 1 倍增到 blockDim.x
//   第 k 轮只有满足 t % stride == 0 的线程活跃
// ---------------------------------------------------------------------------
__global__ void reduceInterleaved2M(const float* __restrict__ input,
                                    float* __restrict__ output,
                                    unsigned int n)
{
    __shared__ float partialSum[2 * BLOCK_SIZE];

    const unsigned int t     = threadIdx.x;
    const unsigned int start = 2 * blockIdx.x * blockDim.x;
    const unsigned int i0    = start + t;
    const unsigned int i1    = start + blockDim.x + t;

    // 两次全局装载：同一 warp 的 32 个 lane 读连续 128 字节 → 完全合并
    // 越界元素写入恒等元 0.0f，使任意 N 都得到正确结果
    partialSum[t]              = (i0 < n) ? input[i0] : 0.0f;
    partialSum[blockDim.x + t] = (i1 < n) ? input[i1] : 0.0f;

    // 循环开头的 __syncthreads()：既保证初始装载完成，也保证每一轮的读
    // 发生在上一轮的写之后。它位于分歧的 if 之外 → 所有线程都会到达。
    for (unsigned int stride = 1; stride <= blockDim.x; stride *= 2) {
        __syncthreads();
        if (t % stride == 0) {
            partialSum[2 * t] += partialSum[2 * t + stride];
        }
    }

    // 循环结束处不需要屏障：最后一轮只有 t == 0 活跃，
    // partialSum[0] 只可能被线程 0 写，线程 0 读自己的写没有跨线程依赖。
    if (t == 0) output[blockIdx.x] = partialSum[0];
}

// ---------------------------------------------------------------------------
// 变体 1B：Kirk/Hwu 经典交错寻址写法（if (tid % (2*s) == 0)）
// ---------------------------------------------------------------------------
__global__ void reduceInterleavedM(const float* __restrict__ input,
                                   float* __restrict__ output,
                                   unsigned int n)
{
    __shared__ float sdata[BLOCK_SIZE];

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

    sdata[tid] = (i < n) ? input[i] : 0.0f;
    __syncthreads();

    for (unsigned int s = 1; s < blockDim.x; s *= 2) {
        if (tid % (2 * s) == 0) {
            sdata[tid] += sdata[tid + s];
        }
        __syncthreads();   // 屏障在分歧的 if 之外 → 安全
    }

    if (tid == 0) output[blockIdx.x] = sdata[0];
}

// ---------------------------------------------------------------------------
// 主机端辅助：打印逐轮发散表（纯解析计算，不依赖 GPU）
// ---------------------------------------------------------------------------
template <typename F>
static void printDivergence(const char* title, int rounds, F perRound)
{
    printf("\n%s\n", title);
    printf("  轮次  活跃线程数  有活干的warp数  有效lane-op  发射的lane槽位  本行lane利用率\n");
    unsigned long long totalUseful = 0, totalSlots = 0;
    for (int r = 1; r <= rounds; ++r) {
        std::pair<unsigned int, unsigned int> pr = perRound(r);
        const unsigned int active = pr.first;
        const unsigned int warps  = pr.second;
        const unsigned long long slots = (unsigned long long)warps * 32ull;
        totalUseful += active;
        totalSlots  += slots;
        printf("  %4d  %10u  %14u  %11u  %15llu  %13.1f%%\n",
               r, active, warps, active, slots,
               100.0 * (double)active / (double)(slots ? slots : 1));
    }
    printf("  合计  有效 lane-op = %llu，发射 lane 槽位 = %llu，整体 lane 有效率 = %.1f%%\n",
           totalUseful, totalSlots, 100.0 * (double)totalUseful / (double)totalSlots);
}

// ---------------------------------------------------------------------------
// 主机端辅助：跑一个 kernel、计时、回传部分和、与 CPU 参考值对比
// ---------------------------------------------------------------------------
typedef void (*KernelFn)(const float*, float*, unsigned int);

static double benchmark(const char* label, KernelFn kfn, unsigned int elemsPerBlock,
                        const float* d_in, unsigned int n, double ref)
{
    const unsigned int numBlocks = (n + elemsPerBlock - 1) / elemsPerBlock;
    float* d_partial = nullptr;
    CUDA_CHECK(cudaMalloc(&d_partial, (size_t)numBlocks * sizeof(float)));

    const int warmup = 3, iters = 20;
    for (int i = 0; i < warmup; ++i) {
        kfn<<<numBlocks, BLOCK_SIZE>>>(d_in, d_partial, n);
    }
    CUDA_CHECK(cudaGetLastError());
    CUDA_CHECK(cudaDeviceSynchronize());

    cudaEvent_t t0, t1;
    CUDA_CHECK(cudaEventCreate(&t0));
    CUDA_CHECK(cudaEventCreate(&t1));
    CUDA_CHECK(cudaEventRecord(t0));
    for (int i = 0; i < iters; ++i) {
        kfn<<<numBlocks, BLOCK_SIZE>>>(d_in, d_partial, n);
    }
    CUDA_CHECK(cudaEventRecord(t1));
    CUDA_CHECK(cudaEventSynchronize(t1));

    float ms = 0.0f;
    CUDA_CHECK(cudaEventElapsedTime(&ms, t0, t1));
    ms /= (float)iters;

    // 收尾：部分和回传主机求和（不计入 kernel 计时）
    float* h_partial = (float*)malloc((size_t)numBlocks * sizeof(float));
    if (!h_partial) { fprintf(stderr, "host malloc failed\n"); exit(EXIT_FAILURE); }
    CUDA_CHECK(cudaMemcpy(h_partial, d_partial, (size_t)numBlocks * sizeof(float),
                          cudaMemcpyDeviceToHost));
    double sum = 0.0;
    for (unsigned int i = 0; i < numBlocks; ++i) sum += (double)h_partial[i];

    const double gbs    = (double)n * (double)sizeof(float) / (ms * 1.0e-3) / 1.0e9;
    const double relErr = (ref != 0.0) ? fabs(sum - ref) / fabs(ref) : fabs(sum - ref);
    printf("%-38s blocks=%-7u 时间=%8.4f ms  有效带宽=%7.1f GB/s  相对误差=%.2e  %s\n",
           label, numBlocks, ms, gbs, relErr, (relErr < 1e-3) ? "PASS" : "FAIL");

    free(h_partial);
    CUDA_CHECK(cudaFree(d_partial));
    CUDA_CHECK(cudaEventDestroy(t0));
    CUDA_CHECK(cudaEventDestroy(t1));
    return ms;
}

int main(int argc, char** argv)
{
    const unsigned int n = (argc > 1)
        ? (unsigned int)strtoul(argv[1], nullptr, 10)
        : 100000000u;
    const size_t bytes = (size_t)n * sizeof(float);

    printf("N = %u 个 float（%.1f MB）\n", n, (double)bytes / (1024.0 * 1024.0));

    cudaDeviceProp prop;
    CUDA_CHECK(cudaGetDeviceProperties(&prop, 0));
    const double theoGBs = 2.0 * (double)prop.memoryClockRate * 1.0e3
                         * ((double)prop.memoryBusWidth / 8.0) / 1.0e9;
    printf("GPU: %s  SM 数 = %d  每 SM 最大线程 = %d  共享内存/SM = %zu B  理论带宽 = %.1f GB/s\n",
           prop.name, prop.multiProcessorCount, prop.maxThreadsPerMultiProcessor,
           prop.sharedMemPerMultiprocessor, theoGBs);
    printf("理论最优时间 = %.4f ms（= N*4 字节 / 理论带宽）\n\n",
           (double)bytes / (theoGBs * 1.0e9) * 1.0e3);

    // ---- 主机端准备数据与 CPU 参考值 ----
    float* h_in = (float*)malloc(bytes);
    if (!h_in) { fprintf(stderr, "host malloc failed\n"); return EXIT_FAILURE; }
    srand(12345);
    double ref = 0.0;
    for (unsigned int i = 0; i < n; ++i) {
        h_in[i] = (float)rand() / (float)RAND_MAX;   // [0,1) 均匀分布
        ref += (double)h_in[i];                      // 用 double 累加作为真值
    }

    float* d_in = nullptr;
    CUDA_CHECK(cudaMalloc(&d_in, bytes));
    CUDA_CHECK(cudaMemcpy(d_in, h_in, bytes, cudaMemcpyHostToDevice));

    // ---- 计时与校验 ----
    const double t1A = benchmark("变体 1A 交错寻址 (2*BLOCK 共享内存)",
                                 reduceInterleaved2M, 2 * BLOCK_SIZE, d_in, n, ref);
    const double t1B = benchmark("变体 1B 交错寻址 (Kirk 形式)",
                                 reduceInterleavedM, BLOCK_SIZE, d_in, n, ref);
    printf("\n变体 1B / 变体 1A 时间比 = %.3f\n", t1B / t1A);

    // ---- 逐轮发散表（解析计算）----
    printDivergence("【变体 1A】逐轮活跃线程与发散度（blockDim = 1024 = 32 个 warp）",
                    11, [](int r) {
        const unsigned int stride = 1u << (r - 1);
        unsigned int active = 0, warps = 0;
        for (unsigned int w = 0; w < BLOCK_SIZE / 32; ++w) {
            bool any = false;
            for (unsigned int l = 0; l < 32; ++l) {
                const unsigned int t = w * 32 + l;
                if (t % stride == 0) { ++active; any = true; }
            }
            if (any) ++warps;
        }
        return std::make_pair(active, warps);
    });
    printDivergence("【变体 1B】逐轮活跃线程与发散度（blockDim = 1024 = 32 个 warp）",
                    10, [](int r) {
        const unsigned int s = 1u << (r - 1);
        unsigned int active = 0, warps = 0;
        for (unsigned int w = 0; w < BLOCK_SIZE / 32; ++w) {
            bool any = false;
            for (unsigned int l = 0; l < 32; ++l) {
                const unsigned int t = w * 32 + l;
                if (t % (2 * s) == 0) { ++active; any = true; }
            }
            if (any) ++warps;
        }
        return std::make_pair(active, warps);
    });

    CUDA_CHECK(cudaFree(d_in));
    free(h_in);
    return EXIT_SUCCESS;
}

// 文件: reduce_v2_sequential.cu
// 编译: nvcc -O3 -arch=sm_80 reduce_v2_sequential.cu -o reduce_v2
// 运行: ./reduce_v2                默认 N = 100000000
//       ./reduce_v2 67108864       自定义 N
//
// 本程序把"交错寻址"与"顺序寻址"放在同一套计时/校验框架里正面比较：
//   reduceInterleaved：交错寻址（示例 1 的变体 1A，作为基线）
//   reduceSequential ：顺序寻址，每线程先合并 2 个元素，
//                      for (s = blockDim.x/2; s > 0; s >>= 1) if (tid < s)

#include <cstdio>
#include <cstdlib>
#include <cmath>
#include <utility>
#include <cuda_runtime.h>

#define BLOCK_SIZE 1024

#define CUDA_CHECK(call)                                                        \
    do {                                                                        \
        cudaError_t err__ = (call);                                             \
        if (err__ != cudaSuccess) {                                             \
            fprintf(stderr, "CUDA error %s:%d: %s\n", __FILE__, __LINE__,       \
                    cudaGetErrorString(err__));                                 \
            exit(EXIT_FAILURE);                                                 \
        }                                                                       \
    } while (0)

// ------------------------------- 基线：交错寻址 -------------------------------
__global__ void reduceInterleaved(const float* __restrict__ input,
                                  float* __restrict__ output,
                                  unsigned int n)
{
    __shared__ float partialSum[2 * BLOCK_SIZE];
    const unsigned int t     = threadIdx.x;
    const unsigned int start = 2 * blockIdx.x * blockDim.x;
    const unsigned int i0    = start + t;
    const unsigned int i1    = start + blockDim.x + t;

    partialSum[t]              = (i0 < n) ? input[i0] : 0.0f;
    partialSum[blockDim.x + t] = (i1 < n) ? input[i1] : 0.0f;

    for (unsigned int stride = 1; stride <= blockDim.x; stride *= 2) {
        __syncthreads();
        if (t % stride == 0) {
            partialSum[2 * t] += partialSum[2 * t + stride];
        }
    }
    if (t == 0) output[blockIdx.x] = partialSum[0];
}

// ------------------------------- 顺序寻址 -------------------------------
// 每线程先从全局内存装载 2 个元素并在寄存器里相加，再进入共享内存树。
// 活跃集合 = 前 s 个连续线程 → warp 内要么全活跃要么全不活跃；地址连续 → 无 bank conflict
__global__ void reduceSequential(const float* __restrict__ input,
                                 float* __restrict__ output,
                                 unsigned int n)
{
    __shared__ float sdata[BLOCK_SIZE];

    const unsigned int tid  = threadIdx.x;
    const unsigned int base = blockIdx.x * (blockDim.x * 2);
    const unsigned int i0   = base + tid;
    const unsigned int i1   = base + blockDim.x + tid;

    const float v0 = (i0 < n) ? input[i0] : 0.0f;
    const float v1 = (i1 < n) ? input[i1] : 0.0f;

    sdata[tid] = v0 + v1;      // 每线程一个元素写共享内存，地址 = tid → 无冲突
    __syncthreads();

    for (unsigned int s = blockDim.x / 2; s > 0; s >>= 1) {
        if (tid < s) {
            sdata[tid] += sdata[tid + s];
        }
        __syncthreads();       // 位于分歧的 if 之外 → 所有线程都会到达
    }

    if (tid == 0) output[blockIdx.x] = sdata[0];
}

// ------------------------------- 主机端工具 -------------------------------
template <typename F>
static void printDivergence(const char* title, int rounds, F perRound)
{
    printf("\n%s\n", title);
    printf("  轮次  活跃线程数  有活干的warp数  有效lane-op  发射的lane槽位  本行lane利用率\n");
    unsigned long long totalUseful = 0, totalSlots = 0;
    for (int r = 1; r <= rounds; ++r) {
        std::pair<unsigned int, unsigned int> pr = perRound(r);
        const unsigned int active = pr.first;
        const unsigned int warps  = pr.second;
        const unsigned long long slots = (unsigned long long)warps * 32ull;
        totalUseful += active;
        totalSlots  += slots;
        printf("  %4d  %10u  %14u  %11u  %15llu  %13.1f%%\n",
               r, active, warps, active, slots,
               100.0 * (double)active / (double)(slots ? slots : 1));
    }
    printf("  合计  有效 lane-op = %llu，发射 lane 槽位 = %llu，整体 lane 有效率 = %.1f%%\n",
           totalUseful, totalSlots, 100.0 * (double)totalUseful / (double)totalSlots);
}

typedef void (*KernelFn)(const float*, float*, unsigned int);

static double benchmark(const char* label, KernelFn kfn, unsigned int elemsPerBlock,
                        const float* d_in, unsigned int n, double ref)
{
    const unsigned int numBlocks = (n + elemsPerBlock - 1) / elemsPerBlock;
    float* d_partial = nullptr;
    CUDA_CHECK(cudaMalloc(&d_partial, (size_t)numBlocks * sizeof(float)));

    for (int i = 0; i < 3; ++i) kfn<<<numBlocks, BLOCK_SIZE>>>(d_in, d_partial, n);
    CUDA_CHECK(cudaGetLastError());
    CUDA_CHECK(cudaDeviceSynchronize());

    cudaEvent_t t0, t1;
    CUDA_CHECK(cudaEventCreate(&t0));
    CUDA_CHECK(cudaEventCreate(&t1));
    CUDA_CHECK(cudaEventRecord(t0));
    for (int i = 0; i < 20; ++i) kfn<<<numBlocks, BLOCK_SIZE>>>(d_in, d_partial, n);
    CUDA_CHECK(cudaEventRecord(t1));
    CUDA_CHECK(cudaEventSynchronize(t1));

    float ms = 0.0f;
    CUDA_CHECK(cudaEventElapsedTime(&ms, t0, t1));
    ms /= 20.0f;

    float* h_partial = (float*)malloc((size_t)numBlocks * sizeof(float));
    if (!h_partial) { fprintf(stderr, "host malloc failed\n"); exit(EXIT_FAILURE); }
    CUDA_CHECK(cudaMemcpy(h_partial, d_partial, (size_t)numBlocks * sizeof(float),
                          cudaMemcpyDeviceToHost));
    double sum = 0.0;
    for (unsigned int i = 0; i < numBlocks; ++i) sum += (double)h_partial[i];

    const double gbs    = (double)n * (double)sizeof(float) / (ms * 1.0e-3) / 1.0e9;
    const double relErr = (ref != 0.0) ? fabs(sum - ref) / fabs(ref) : fabs(sum - ref);
    printf("%-40s blocks=%-7u 时间=%8.4f ms  有效带宽=%7.1f GB/s  相对误差=%.2e  %s\n",
           label, numBlocks, ms, gbs, relErr, (relErr < 1e-3) ? "PASS" : "FAIL");

    free(h_partial);
    CUDA_CHECK(cudaFree(d_partial));
    CUDA_CHECK(cudaEventDestroy(t0));
    CUDA_CHECK(cudaEventDestroy(t1));
    return ms;
}

int main(int argc, char** argv)
{
    const unsigned int n = (argc > 1)
        ? (unsigned int)strtoul(argv[1], nullptr, 10)
        : 100000000u;
    const size_t bytes = (size_t)n * sizeof(float);

    cudaDeviceProp prop;
    CUDA_CHECK(cudaGetDeviceProperties(&prop, 0));
    const double theoGBs = 2.0 * (double)prop.memoryClockRate * 1.0e3
                         * ((double)prop.memoryBusWidth / 8.0) / 1.0e9;
    printf("GPU: %s  理论带宽 = %.1f GB/s  N = %u（%.1f MB）  理论下限 = %.4f ms\n\n",
           prop.name, theoGBs, n, (double)bytes / (1024.0 * 1024.0),
           (double)bytes / (theoGBs * 1.0e9) * 1.0e3);

    float* h_in = (float*)malloc(bytes);
    if (!h_in) { fprintf(stderr, "host malloc failed\n"); return EXIT_FAILURE; }
    srand(2025);
    double ref = 0.0;
    for (unsigned int i = 0; i < n; ++i) {
        h_in[i] = (float)rand() / (float)RAND_MAX;
        ref += (double)h_in[i];
    }

    float* d_in = nullptr;
    CUDA_CHECK(cudaMalloc(&d_in, bytes));
    CUDA_CHECK(cudaMemcpy(d_in, h_in, bytes, cudaMemcpyHostToDevice));

    const double tV1 = benchmark("V1 交错寻址   (t % stride == 0)",
                                 reduceInterleaved, 2 * BLOCK_SIZE, d_in, n, ref);
    const double tV2 = benchmark("V2 顺序寻址   (tid < s)",
                                 reduceSequential, 2 * BLOCK_SIZE, d_in, n, ref);
    printf("\n顺序寻址 / 交错寻址 时间比 = %.3f（小于 1 表示顺序寻址更快）\n", tV2 / tV1);
    printf("带宽提升 = %.1f GB/s → %.1f GB/s\n",
           (double)n * 4.0 / (tV1 * 1.0e-3) / 1.0e9,
           (double)n * 4.0 / (tV2 * 1.0e-3) / 1.0e9);

    printDivergence("【顺序寻址】逐轮活跃线程与发散度（blockDim = 1024 = 32 个 warp）",
                    10, [](int r) {
        const unsigned int s = BLOCK_SIZE >> r;   // r=1 → 512, r=10 → 1
        unsigned int active = 0, warps = 0;
        for (unsigned int w = 0; w < BLOCK_SIZE / 32; ++w) {
            bool any = false;
            for (unsigned int l = 0; l < 32; ++l) {
                const unsigned int t = w * 32 + l;
                if (t < s) { ++active; any = true; }
            }
            if (any) ++warps;
        }
        return std::make_pair(active, warps);
    });

    CUDA_CHECK(cudaFree(d_in));
    free(h_in);
    return EXIT_SUCCESS;
}

// 文件: reduce_v3_shfl.cu
// 编译: nvcc -O3 -arch=sm_80 reduce_v3_shfl.cu -o reduce_v3
// 运行: ./reduce_v3              默认 N = 100000000，同时测试 BLOCK = 256 与 512
//       ./reduce_v3 33554432     自定义 N
//
// 本程序在同一个进程里依次测量四个 kernel，直接给出对比表格：
//   V1 交错寻址（示例 1 变体 1A）
//   V2 顺序寻址（示例 2）
//   V3 粗化 + 5 轮 __shfl_down_sync + 每个 warp 只写 1 个共享内存字 + 每 block 一次 atomicAdd
//      BLOCK = 256 与 BLOCK = 512 各测一遍，并用 cudaOccupancyMaxActiveBlocksPerMultiprocessor
//      计算理论占用率、用 cudaFuncGetAttributes 打印寄存器用量。

#include <cstdio>
#include <cstdlib>
#include <cmath>
#include <cuda_runtime.h>

#define BLOCK_SIZE 1024          // V1 / V2 使用

#define CUDA_CHECK(call)                                                        \
    do {                                                                        \
        cudaError_t err__ = (call);                                             \
        if (err__ != cudaSuccess) {                                             \
            fprintf(stderr, "CUDA error %s:%d: %s\n", __FILE__, __LINE__,       \
                    cudaGetErrorString(err__));                                 \
            exit(EXIT_FAILURE);                                                 \
        }                                                                       \
    } while (0)

// ------------------------------- V1：交错寻址（基线） -------------------------------
__global__ void reduceInterleaved(const float* __restrict__ input,
                                  float* __restrict__ output, unsigned int n)
{
    __shared__ float partialSum[2 * BLOCK_SIZE];
    const unsigned int t     = threadIdx.x;
    const unsigned int start = 2 * blockIdx.x * blockDim.x;
    const unsigned int i0    = start + t;
    const unsigned int i1    = start + blockDim.x + t;

    partialSum[t]              = (i0 < n) ? input[i0] : 0.0f;
    partialSum[blockDim.x + t] = (i1 < n) ? input[i1] : 0.0f;

    for (unsigned int stride = 1; stride <= blockDim.x; stride *= 2) {
        __syncthreads();
        if (t % stride == 0) partialSum[2 * t] += partialSum[2 * t + stride];
    }
    if (t == 0) output[blockIdx.x] = partialSum[0];
}

// ------------------------------- V2：顺序寻址（基线） -------------------------------
__global__ void reduceSequential(const float* __restrict__ input,
                                 float* __restrict__ output, unsigned int n)
{
    __shared__ float sdata[BLOCK_SIZE];
    const unsigned int tid  = threadIdx.x;
    const unsigned int base = blockIdx.x * (blockDim.x * 2);
    const unsigned int i0   = base + tid;
    const unsigned int i1   = base + blockDim.x + tid;

    const float v0 = (i0 < n) ? input[i0] : 0.0f;
    const float v1 = (i1 < n) ? input[i1] : 0.0f;

    sdata[tid] = v0 + v1;
    __syncthreads();

    for (unsigned int s = blockDim.x / 2; s > 0; s >>= 1) {
        if (tid < s) sdata[tid] += sdata[tid + s];
        __syncthreads();
    }
    if (tid == 0) output[blockIdx.x] = sdata[0];
}

// ------------------------------- V3：粗化 + shuffle + 原子收尾 -------------------------------
// BLOCK 为编译期常量（模板参数），便于让编译器把 warpSums 的大小与 5 轮 shuffle 完全展开。
template <unsigned int BLOCK>
__global__ void reduceCoarsenedShfl(const float* __restrict__ input,
                                    float* __restrict__ output,   // 单个标量，启动前须清零
                                    unsigned int n)
{
    __shared__ float warpSums[BLOCK / 32];

    const unsigned int tid  = threadIdx.x;
    const unsigned int lane = tid & 31u;          // warp 内 lane 编号 0..31
    const unsigned int wid  = tid >> 5;           // block 内 warp 编号
    const unsigned int G    = gridDim.x * BLOCK;  // 网格总线程数 = 网格-步长循环的步长
    const unsigned int base = blockIdx.x * BLOCK + tid;

    // ---- ① 线程粗化：每个线程用网格-步长循环累加数百个元素 ----
    // 4 个独立累加器把 FADD 依赖链长度除以 4，提高指令级并行（ILP）
    float s0 = 0.0f, s1 = 0.0f, s2 = 0.0f, s3 = 0.0f;
    unsigned int i = base;
    for (; i + 3u * G < n; i += 4u * G) {
        s0 += input[i];
        s1 += input[i + G];
        s2 += input[i + 2u * G];
        s3 += input[i + 3u * G];
    }
    for (; i < n; i += G) s0 += input[i];
    float sum = (s0 + s1) + (s2 + s3);

    // ---- ② warp 内归约：5 轮 __shfl_down_sync，不需要 __syncthreads() ----
    #pragma unroll
    for (int offset = 16; offset > 0; offset >>= 1) {
        sum += __shfl_down_sync(0xffffffffu, sum, offset);
    }

    // ---- ③ 每个 warp 的 lane 0 写 1 个 float 到共享内存 ----
    // 一个 warp 只有一条活跃 lane 执行写操作 → 不可能发生 bank conflict
    if (lane == 0) warpSums[wid] = sum;
    __syncthreads();        // 唯一的屏障：保证所有 warp 的部分和都已落盘

    // ---- ④ block 内最后一个 warp（warp 0）完成最终归约，同样用 shuffle ----
    if (wid == 0) {
        const unsigned int nwarps = BLOCK / 32u;
        float v = (lane < nwarps) ? warpSums[lane] : 0.0f;   // 越界 lane 用恒等元 0.0f
        #pragma unroll
        for (int offset = 16; offset > 0; offset >>= 1) {
            v += __shfl_down_sync(0xffffffffu, v, offset);
        }
        if (lane == 0) atomicAdd(output, v);   // 每个 block 只发一条原子操作
    }
}

// ------------------------------- 主机端工具 -------------------------------
static double cpuReference(const float* h, unsigned int n)
{
    double s = 0.0;
    for (unsigned int i = 0; i < n; ++i) s += (double)h[i];
    return s;
}

// 通用基准：kernel 输出"每 block 一个部分和"，本函数负责计时 + 主机收尾 + 校验
typedef void (*KernelFn)(const float*, float*, unsigned int);

static double benchPartitioned(const char* label, KernelFn kfn, unsigned int block,
                               unsigned int elemsPerBlock, const float* d_in,
                               unsigned int n, double ref)
{
    const unsigned int numBlocks = (n + elemsPerBlock - 1) / elemsPerBlock;
    float* d_partial = nullptr;
    CUDA_CHECK(cudaMalloc(&d_partial, (size_t)numBlocks * sizeof(float)));

    for (int i = 0; i < 3; ++i) kfn<<<numBlocks, block>>>(d_in, d_partial, n);
    CUDA_CHECK(cudaGetLastError());
    CUDA_CHECK(cudaDeviceSynchronize());

    cudaEvent_t t0, t1;
    CUDA_CHECK(cudaEventCreate(&t0));
    CUDA_CHECK(cudaEventCreate(&t1));
    CUDA_CHECK(cudaEventRecord(t0));
    for (int i = 0; i < 20; ++i) kfn<<<numBlocks, block>>>(d_in, d_partial, n);
    CUDA_CHECK(cudaEventRecord(t1));
    CUDA_CHECK(cudaEventSynchronize(t1));

    float ms = 0.0f;
    CUDA_CHECK(cudaEventElapsedTime(&ms, t0, t1));
    ms /= 20.0f;

    float* h_partial = (float*)malloc((size_t)numBlocks * sizeof(float));
    CUDA_CHECK(cudaMemcpy(h_partial, d_partial, (size_t)numBlocks * sizeof(float),
                          cudaMemcpyDeviceToHost));
    double sum = 0.0;
    for (unsigned int i = 0; i < numBlocks; ++i) sum += (double)h_partial[i];

    const double gbs    = (double)n * 4.0 / (ms * 1.0e-3) / 1.0e9;
    const double relErr = fabs(sum - ref) / fabs(ref);
    printf("%-34s blocks=%-7u 时间=%8.4f ms  带宽=%7.1f GB/s  误差=%.2e  %s\n",
           label, numBlocks, ms, gbs, relErr, (relErr < 1e-3) ? "PASS" : "FAIL");

    free(h_partial);
    CUDA_CHECK(cudaFree(d_partial));
    CUDA_CHECK(cudaEventDestroy(t0));
    CUDA_CHECK(cudaEventDestroy(t1));
    return ms;
}

// V3 专用基准：输出是单个标量（atomicAdd 收尾），需要每轮清零
template <unsigned int BLOCK>
static double benchAtomic(const char* label, const float* d_in, unsigned int n,
                          double ref, int gridBlocks)
{
    float* d_out = nullptr;
    CUDA_CHECK(cudaMalloc(&d_out, sizeof(float)));

    // 计时：连续 20 次启动，不在计时区间里做 memset（累加值不影响带宽测量）
    for (int i = 0; i < 3; ++i) reduceCoarsenedShfl<BLOCK><<<gridBlocks, BLOCK>>>(d_in, d_out, n);
    CUDA_CHECK(cudaGetLastError());
    CUDA_CHECK(cudaDeviceSynchronize());

    cudaEvent_t t0, t1;
    CUDA_CHECK(cudaEventCreate(&t0));
    CUDA_CHECK(cudaEventCreate(&t1));
    CUDA_CHECK(cudaEventRecord(t0));
    for (int i = 0; i < 20; ++i) reduceCoarsenedShfl<BLOCK><<<gridBlocks, BLOCK>>>(d_in, d_out, n);
    CUDA_CHECK(cudaEventRecord(t1));
    CUDA_CHECK(cudaEventSynchronize(t1));
    float ms = 0.0f;
    CUDA_CHECK(cudaEventElapsedTime(&ms, t0, t1));
    ms /= 20.0f;

    // 校验：清零后单独跑一次
    CUDA_CHECK(cudaMemset(d_out, 0, sizeof(float)));
    reduceCoarsenedShfl<BLOCK><<<gridBlocks, BLOCK>>>(d_in, d_out, n);
    CUDA_CHECK(cudaGetLastError());
    float h_out = 0.0f;
    CUDA_CHECK(cudaMemcpy(&h_out, d_out, sizeof(float), cudaMemcpyDeviceToHost));

    const double gbs    = (double)n * 4.0 / (ms * 1.0e-3) / 1.0e9;
    const double relErr = fabs((double)h_out - ref) / fabs(ref);
    printf("%-34s blocks=%-7d 时间=%8.4f ms  带宽=%7.1f GB/s  误差=%.2e  %s\n",
           label, gridBlocks, ms, gbs, relErr, (relErr < 1e-3) ? "PASS" : "FAIL");

    CUDA_CHECK(cudaFree(d_out));
    CUDA_CHECK(cudaEventDestroy(t0));
    CUDA_CHECK(cudaEventDestroy(t1));
    return ms;
}

// 打印占用率与寄存器信息
template <unsigned int BLOCK>
static void reportOccupancy(const char* label, const cudaDeviceProp& prop)
{
    int blocksPerSm = 0;
    CUDA_CHECK(cudaOccupancyMaxActiveBlocksPerMultiprocessor(
        &blocksPerSm, reduceCoarsenedShfl<BLOCK>, BLOCK, 0));

    cudaFuncAttributes attr;
    CUDA_CHECK(cudaFuncGetAttributes(&attr, reduceCoarsenedShfl<BLOCK>));

    const int threadsPerSm = blocksPerSm * (int)BLOCK;
    const double occ = 100.0 * (double)threadsPerSm / (double)prop.maxThreadsPerMultiProcessor;
    printf("  %s: 寄存器/线程 = %d，静态共享内存 = %zu B，"
           "每 SM 驻留 block = %d，每 SM 线程 = %d / %d，占用率 = %.1f%%\n",
           label, attr.numRegs, attr.sharedSizeBytes, blocksPerSm,
           threadsPerSm, prop.maxThreadsPerMultiProcessor, occ);
}

int main(int argc, char** argv)
{
    const unsigned int n = (argc > 1)
        ? (unsigned int)strtoul(argv[1], nullptr, 10)
        : 100000000u;
    const size_t bytes = (size_t)n * sizeof(float);

    cudaDeviceProp prop;
    CUDA_CHECK(cudaGetDeviceProperties(&prop, 0));
    const double theoGBs = 2.0 * (double)prop.memoryClockRate * 1.0e3
                         * ((double)prop.memoryBusWidth / 8.0) / 1.0e9;
    printf("GPU: %s  理论带宽 = %.1f GB/s  N = %u（%.1f MB）  理论下限 = %.4f ms\n\n",
           prop.name, theoGBs, n, (double)bytes / (1024.0 * 1024.0),
           (double)bytes / (theoGBs * 1.0e9) * 1.0e3);

    float* h_in = (float*)malloc(bytes);
    srand(4096);
    for (unsigned int i = 0; i < n; ++i) h_in[i] = (float)rand() / (float)RAND_MAX;
    const double ref = cpuReference(h_in, n);

    float* d_in = nullptr;
    CUDA_CHECK(cudaMalloc(&d_in, bytes));
    CUDA_CHECK(cudaMemcpy(d_in, h_in, bytes, cudaMemcpyHostToDevice));

    printf("================ 基线：V1 / V2 ================\n");
    const double tV1 = benchPartitioned("V1 交错寻址", reduceInterleaved,
                                        BLOCK_SIZE, 2 * BLOCK_SIZE, d_in, n, ref);
    const double tV2 = benchPartitioned("V2 顺序寻址", reduceSequential,
                                        BLOCK_SIZE, 2 * BLOCK_SIZE, d_in, n, ref);

    printf("\n================ V3：粗化 + shuffle + 原子 ================\n");
    reportOccupancy<256>("BLOCK=256", prop);
    reportOccupancy<512>("BLOCK=512", prop);

    // 网格规模 = 每 SM 可驻留的 block 数 × SM 数（整数倍满载，令粗化因子最大化）
    int bps256 = 0, bps512 = 0;
    CUDA_CHECK(cudaOccupancyMaxActiveBlocksPerMultiprocessor(
        &bps256, reduceCoarsenedShfl<256>, 256, 0));
    CUDA_CHECK(cudaOccupancyMaxActiveBlocksPerMultiprocessor(
        &bps512, reduceCoarsenedShfl<512>, 512, 0));
    const int grid256 = bps256 * prop.multiProcessorCount;
    const int grid512 = bps512 * prop.multiProcessorCount;
    printf("  网格规模: BLOCK=256 → %d 个 block（每线程约 %.0f 个元素）；"
           "BLOCK=512 → %d 个 block（每线程约 %.0f 个元素）\n",
           grid256, (double)n / ((double)grid256 * 256.0),
           grid512, (double)n / ((double)grid512 * 512.0));

    const double t256 = benchAtomic<256>("V3 BLOCK=256 粗化+shuffle", d_in, n, ref, grid256);
    const double t512 = benchAtomic<512>("V3 BLOCK=512 粗化+shuffle", d_in, n, ref, grid512);

    printf("\n================ 汇总 ================\n");
    printf("  理论下限                    : %8.4f ms  (%.1f GB/s, 100%%)\n",
           (double)bytes / (theoGBs * 1.0e9) * 1.0e3, theoGBs);
    printf("  V1 交错寻址                 : %8.4f ms  (%.1f GB/s, %.1f%%)  相对 V1 = %.2fx\n",
           tV1, (double)n * 4.0 / (tV1 * 1.0e-3) / 1.0e9,
           100.0 * (double)n * 4.0 / (tV1 * 1.0e-3) / (theoGBs * 1.0e9), 1.0);
    printf("  V2 顺序寻址                 : %8.4f ms  (%.1f GB/s, %.1f%%)  相对 V1 = %.2fx\n",
           tV2, (double)n * 4.0 / (tV2 * 1.0e-3) / 1.0e9,
           100.0 * (double)n * 4.0 / (tV2 * 1.0e-3) / (theoGBs * 1.0e9), tV1 / tV2);
    printf("  V3 BLOCK=256                : %8.4f ms  (%.1f GB/s, %.1f%%)  相对 V1 = %.2fx\n",
           t256, (double)n * 4.0 / (t256 * 1.0e-3) / 1.0e9,
           100.0 * (double)n * 4.0 / (t256 * 1.0e-3) / (theoGBs * 1.0e9), tV1 / t256);
    printf("  V3 BLOCK=512                : %8.4f ms  (%.1f GB/s, %.1f%%)  相对 V1 = %.2fx\n",
           t512, (double)n * 4.0 / (t512 * 1.0e-3) / 1.0e9,
           100.0 * (double)n * 4.0 / (t512 * 1.0e-3) / (theoGBs * 1.0e9), tV1 / t512);

    CUDA_CHECK(cudaFree(d_in));
    free(h_in);
    return EXIT_SUCCESS;
}

// 文件: reduce_v4_hierarchical.cu
// 编译: nvcc -O3 -arch=sm_80 reduce_v4_hierarchical.cu -o reduce_v4
// 运行: ./reduce_v4                   默认 N = 100000003（故意不是 4 的倍数，考察边界处理）
//       ./reduce_v4 100000000         自定义 N
//       ./reduce_v4 100000000 stress  使用量级悬殊的数据（一半 1e3、一半 1e-3）
//
// 本程序演示两种"grid 级收尾"方式，并对比它们的正确性、性能与可复现性：
//   (a) 多阶段 kernel：阶段 1 产出 864 个部分和，阶段 2 产出 4 个，阶段 3 产出 1 个
//       —— 归约树形状固定，结果逐位可复现
//   (b) 单阶段 kernel + atomicAdd：每个 block 一条原子操作
//       —— 浮点加法不满足结合律，更新顺序不定 → 结果可能不可复现
// 阶段 1 用 float4 向量装载 + 标量尾巴处理任意 N（含 n % 4 != 0 的情形）。

#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <cmath>
#include <vector>
#include <cuda_runtime.h>

#define CUDA_CHECK(call)                                                        \
    do {                                                                        \
        cudaError_t err__ = (call);                                             \
        if (err__ != cudaSuccess) {                                             \
            fprintf(stderr, "CUDA error %s:%d: %s\n", __FILE__, __LINE__,       \
                    cudaGetErrorString(err__));                                 \
            exit(EXIT_FAILURE);                                                 \
        }                                                                       \
    } while (0)

// 前置声明：下面 4(c) 用这个"每线程一条 atomicAdd"的反例 kernel 做对照实验，
// 它的定义见本节末尾的"反例代码"块（放在 main 之后只是为了先讲清楚主线）。
__global__ void reduceEveryThreadAtomic(const float* __restrict__ in,
                                        unsigned int n, float* __restrict__ out);

// ---------------------------------------------------------------------------
// 通用的"分段归约"kernel：
//   USE_ATOMIC == false → 每个 block 写 out[blockIdx.x]（多阶段用，确定性）
//   USE_ATOMIC == true  → 每个 block 累加到 out[0]（单阶段用，需提前清零）
// 主体用 float4 向量装载并对齐到 16 字节；n % 4 的余数由 block 0 用标量循环处理。
// ---------------------------------------------------------------------------
template <unsigned int BLOCK, bool USE_ATOMIC>
__global__ void reduceStage(const float* __restrict__ in, unsigned int n,
                            float* __restrict__ out)
{
    __shared__ float warpSums[BLOCK / 32];

    const unsigned int tid  = threadIdx.x;
    const unsigned int lane = tid & 31u;
    const unsigned int wid  = tid >> 5;
    const unsigned int G    = gridDim.x * BLOCK;
    const unsigned int n4   = n >> 2;                    // float4 的个数
    const float4* __restrict__ in4 = reinterpret_cast<const float4*>(in);

    // ---- ① 向量化 + 粗化：主循环遍历完整的 float4 ----
    float sum = 0.0f;
    for (unsigned int i = blockIdx.x * BLOCK + tid; i < n4; i += G) {
        const float4 v = in4[i];
        sum += (v.x + v.y) + (v.z + v.w);
    }

    // ---- ② 标量尾巴：n % 4 个元素由 block 0 处理，避免重复计数 ----
    if (blockIdx.x == 0) {
        for (unsigned int i = (n4 << 2) + tid; i < n; i += BLOCK) {
            sum += in[i];
        }
    }

    // ---- ③ warp 内 shuffle 归约 ----
    #pragma unroll
    for (int offset = 16; offset > 0; offset >>= 1) {
        sum += __shfl_down_sync(0xffffffffu, sum, offset);
    }
    if (lane == 0) warpSums[wid] = sum;
    __syncthreads();

    // ---- ④ 最后一个 warp 完成 block 内最终归约 ----
    if (wid == 0) {
        const unsigned int nwarps = BLOCK / 32u;
        float v = (lane < nwarps) ? warpSums[lane] : 0.0f;
        #pragma unroll
        for (int offset = 16; offset > 0; offset >>= 1) {
            v += __shfl_down_sync(0xffffffffu, v, offset);
        }
        if (lane == 0) {
            if (USE_ATOMIC) atomicAdd(out, v);
            else            out[blockIdx.x] = v;
        }
    }
}

int main(int argc, char** argv)
{
    const unsigned int n = (argc > 1)
        ? (unsigned int)strtoul(argv[1], nullptr, 10)
        : 100000003u;
    const bool stress = (argc > 2) && (strcmp(argv[2], "stress") == 0);
    const size_t bytes = (size_t)n * sizeof(float);

    cudaDeviceProp prop;
    CUDA_CHECK(cudaGetDeviceProperties(&prop, 0));
    const double theoGBs = 2.0 * (double)prop.memoryClockRate * 1.0e3
                         * ((double)prop.memoryBusWidth / 8.0) / 1.0e9;
    const double tMin = (double)bytes / (theoGBs * 1.0e9) * 1.0e3;
    printf("GPU: %s  理论带宽 = %.1f GB/s  N = %u（%.1f MB）  理论下限 = %.4f ms%s\n\n",
           prop.name, theoGBs, n, (double)bytes / (1024.0 * 1024.0), tMin,
           stress ? "  [stress 数据：量级相差 10^6 倍]" : "");

    // ---- 主机端数据：stress 模式故意让量级悬殊，逼迫浮点顺序误差显形 ----
    float* h_in = (float*)malloc(bytes);
    if (!h_in) { fprintf(stderr, "host malloc failed\n"); return EXIT_FAILURE; }
    srand(31337);
    for (unsigned int i = 0; i < n; ++i) {
        if (stress) h_in[i] = ((i & 1u) == 0u) ? 1.0e3f : 1.0e-3f;
        else        h_in[i] = (float)rand() / (float)RAND_MAX;
    }
    double ref = 0.0;
    for (unsigned int i = 0; i < n; ++i) ref += (double)h_in[i];

    float* d_in = nullptr;
    CUDA_CHECK(cudaMalloc(&d_in, bytes));
    CUDA_CHECK(cudaMemcpy(d_in, h_in, bytes, cudaMemcpyHostToDevice));

    const unsigned int BLOCK = 256;
    int bps = 0;
    CUDA_CHECK(cudaOccupancyMaxActiveBlocksPerMultiprocessor(
        &bps, reduceStage<BLOCK, false>, BLOCK, 0));
    const int maxGrid = bps * prop.multiProcessorCount;
    printf("BLOCK = %u，每 SM 驻留 %d 个 block，最大网格 = %d 个 block（%d 个线程）\n",
           BLOCK, bps, maxGrid, maxGrid * (int)BLOCK);

    // ---- 缓冲区：多阶段版本需要一个 ping-pong 缓冲，长度取 maxGrid ----
    float* d_bufA = nullptr;
    float* d_bufB = nullptr;
    CUDA_CHECK(cudaMalloc(&d_bufA, (size_t)maxGrid * sizeof(float)));
    CUDA_CHECK(cudaMalloc(&d_bufB, (size_t)maxGrid * sizeof(float)));

    cudaEvent_t t0, t1;
    CUDA_CHECK(cudaEventCreate(&t0));
    CUDA_CHECK(cudaEventCreate(&t1));

    // ============================ (a) 多阶段 kernel ============================
    unsigned int count = n;
    const float* src = d_in;
    float* dst = d_bufA;
    int stages = 0;

    // 预热
    {
        unsigned int c = n;
        const float* s = d_in;
        float* d = d_bufA;
        while (c > 1) {
            unsigned int blocks = (c + BLOCK - 1u) / BLOCK;
            if (blocks > (unsigned int)maxGrid) blocks = (unsigned int)maxGrid;
            reduceStage<BLOCK, false><<<blocks, BLOCK>>>(s, c, d);
            c = blocks;
            s = d;
            d = (d == d_bufA) ? d_bufB : d_bufA;
        }
    }
    CUDA_CHECK(cudaGetLastError());
    CUDA_CHECK(cudaDeviceSynchronize());

    CUDA_CHECK(cudaEventRecord(t0));
    for (int rep = 0; rep < 20; ++rep) {
        count = n;
        src = d_in;
        dst = d_bufA;
        stages = 0;
        while (count > 1) {
            unsigned int blocks = (count + BLOCK - 1u) / BLOCK;
            if (blocks > (unsigned int)maxGrid) blocks = (unsigned int)maxGrid;
            reduceStage<BLOCK, false><<<blocks, BLOCK>>>(src, count, dst);
            count = blocks;
            src = dst;
            dst = (dst == d_bufA) ? d_bufB : d_bufA;
            ++stages;
        }
    }
    CUDA_CHECK(cudaEventRecord(t1));
    CUDA_CHECK(cudaEventSynchronize(t1));
    float msMulti = 0.0f;
    CUDA_CHECK(cudaEventElapsedTime(&msMulti, t0, t1));
    msMulti /= 20.0f;

    float h_multi = 0.0f;
    CUDA_CHECK(cudaMemcpy(&h_multi, src, sizeof(float), cudaMemcpyDeviceToHost));
    const double errMulti = fabs((double)h_multi - ref) / fabs(ref);

    // 可复现性：重复 10 次，比较结果的二进制位
    bool deterministic = true;
    float first = h_multi;
    for (int rep = 0; rep < 10; ++rep) {
        unsigned int c = n;
        const float* s = d_in;
        float* d = d_bufA;
        while (c > 1) {
            unsigned int blocks = (c + BLOCK - 1u) / BLOCK;
            if (blocks > (unsigned int)maxGrid) blocks = (unsigned int)maxGrid;
            reduceStage<BLOCK, false><<<blocks, BLOCK>>>(s, c, d);
            c = blocks;
            s = d;
            d = (d == d_bufA) ? d_bufB : d_bufA;
        }
        float h = 0.0f;
        CUDA_CHECK(cudaMemcpy(&h, s, sizeof(float), cudaMemcpyDeviceToHost));
        if (memcmp(&h, &first, sizeof(float)) != 0) deterministic = false;
    }

    printf("(a) 多阶段 kernel（%d 个阶段）  : 时间=%8.4f ms  带宽=%7.1f GB/s  误差=%.2e  %s  逐位可复现=%s\n",
           stages, msMulti, (double)n * 4.0 / (msMulti * 1.0e-3) / 1.0e9, errMulti,
           (errMulti < 1e-3) ? "PASS" : "FAIL", deterministic ? "是" : "否");

    // ============================ (b) 单阶段 + atomicAdd ============================
    float* d_sum = nullptr;
    CUDA_CHECK(cudaMalloc(&d_sum, sizeof(float)));

    for (int rep = 0; rep < 3; ++rep) {
        CUDA_CHECK(cudaMemset(d_sum, 0, sizeof(float)));
        reduceStage<BLOCK, true><<<maxGrid, BLOCK>>>(d_in, n, d_sum);
    }
    CUDA_CHECK(cudaGetLastError());
    CUDA_CHECK(cudaDeviceSynchronize());

    CUDA_CHECK(cudaEventRecord(t0));
    for (int rep = 0; rep < 20; ++rep) {
        reduceStage<BLOCK, true><<<maxGrid, BLOCK>>>(d_in, n, d_sum);
    }
    CUDA_CHECK(cudaEventRecord(t1));
    CUDA_CHECK(cudaEventSynchronize(t1));
    float msAtomic = 0.0f;
    CUDA_CHECK(cudaEventElapsedTime(&msAtomic, t0, t1));
    msAtomic /= 20.0f;

    CUDA_CHECK(cudaMemset(d_sum, 0, sizeof(float)));
    reduceStage<BLOCK, true><<<maxGrid, BLOCK>>>(d_in, n, d_sum);
    float h_atomic = 0.0f;
    CUDA_CHECK(cudaMemcpy(&h_atomic, d_sum, sizeof(float), cudaMemcpyDeviceToHost));
    const double errAtomic = fabs((double)h_atomic - ref) / fabs(ref);

    // 可复现性：跑 20 次，统计有多少次与第一次的二进制位不同
    int different = 0;
    const float firstAtomic = h_atomic;
    for (int rep = 0; rep < 20; ++rep) {
        CUDA_CHECK(cudaMemset(d_sum, 0, sizeof(float)));
        reduceStage<BLOCK, true><<<maxGrid, BLOCK>>>(d_in, n, d_sum);
        float h = 0.0f;
        CUDA_CHECK(cudaMemcpy(&h, d_sum, sizeof(float), cudaMemcpyDeviceToHost));
        if (memcmp(&h, &firstAtomic, sizeof(float)) != 0) ++different;
    }

    printf("(b) 单阶段 + atomicAdd        : 时间=%8.4f ms  带宽=%7.1f GB/s  误差=%.2e  %s  逐位可复现=%s（20 次里 %d 次不同）\n",
           msAtomic, (double)n * 4.0 / (msAtomic * 1.0e-3) / 1.0e9, errAtomic,
           (errAtomic < 1e-3) ? "PASS" : "FAIL",
           (different == 0) ? "是" : "否", different);

    printf("\nCPU 参考值 = %.10e   多阶段结果 = %.10e   原子结果 = %.10e\n",
           ref, (double)h_multi, (double)h_atomic);
    printf("理论下限 = %.4f ms；多阶段达到峰值带宽的 %.1f%%，原子版达到 %.1f%%\n",
           tMin, 100.0 * (double)n * 4.0 / (msMulti * 1.0e-3) / (theoGBs * 1.0e9),
           100.0 * (double)n * 4.0 / (msAtomic * 1.0e-3) / (theoGBs * 1.0e9));

    // ============================ (c) 故意慢的写法：每线程一条 atomicAdd ============================
    // 只在小规模 N 上演示，否则要在 L2 里排几分钟的队
    if (n <= 4000000u) {
        const unsigned int smallN = n;
        float* d_slow = nullptr;
        CUDA_CHECK(cudaMalloc(&d_slow, sizeof(float)));
        CUDA_CHECK(cudaMemset(d_slow, 0, sizeof(float)));
        CUDA_CHECK(cudaEventRecord(t0));
        // 启动 smallN 个线程，每个线程读一个元素后立刻 atomicAdd 到同一个地址
        reduceEveryThreadAtomic<<<(smallN + 255u) / 256u, 256>>>(d_in, smallN, d_slow);
        CUDA_CHECK(cudaEventRecord(t1));
        CUDA_CHECK(cudaEventSynchronize(t1));
        float msSlow = 0.0f;
        CUDA_CHECK(cudaEventElapsedTime(&msSlow, t0, t1));
        float h_slow = 0.0f;
        CUDA_CHECK(cudaMemcpy(&h_slow, d_slow, sizeof(float), cudaMemcpyDeviceToHost));
        printf("(c) 每线程一条 atomicAdd（N=%u）: 时间=%8.4f ms  带宽=%7.1f GB/s  误差=%.2e\n",
               smallN, msSlow, (double)smallN * 4.0 / (msSlow * 1.0e-3) / 1.0e9,
               fabs((double)h_slow - ref) / fabs(ref));
        CUDA_CHECK(cudaFree(d_slow));
    } else {
        printf("(c) 每线程一条 atomicAdd：N 太大（%u > 4000000），按 3e8 次/秒的同址原子吞吐估算，"
               "需要约 %.1f ms，是理论下限的 %.0f 倍，故跳过实测\n",
               n, (double)n / 3.0e8 * 1.0e3, (double)n / 3.0e8 * 1.0e3 / tMin);
    }

    CUDA_CHECK(cudaFree(d_bufA));
    CUDA_CHECK(cudaFree(d_bufB));
    CUDA_CHECK(cudaFree(d_sum));
    CUDA_CHECK(cudaFree(d_in));
    CUDA_CHECK(cudaEventDestroy(t0));
    CUDA_CHECK(cudaEventDestroy(t1));
    free(h_in);
    return EXIT_SUCCESS;
}

// 反例：每个线程读完一个元素就 atomicAdd 到同一个地址。
// 10^8 个线程全部竞争同一个 L2 地址 → 原子操作被串行化，吞吐量约 2-5 × 10^8 次/秒。
__global__ void reduceEveryThreadAtomic(const float* __restrict__ in,
                                        unsigned int n, float* __restrict__ out)
{
    const unsigned int i = blockIdx.x * blockDim.x + threadIdx.x;
    if (i < n) atomicAdd(out, in[i]);
}