Lecture 12: 最终项目 —— 问题定义、方案设计、性能优化与报告 (对应 Final Project)

// 文件: project_template.cu
// 编译: nvcc -O3 -arch=sm_80 -lineinfo project_template.cu -o project_template
// 运行: ./project_template -m 1000 -n 1000 -k 1000 -tile 16 -bs 256 -iters 20
//       ./project_template -m 1 -n 1 -k 1                # 边界: 单个元素
//       ./project_template -m 1025 -n 1023 -k 1000       # 边界: 非 2 的幂 + 非整数倍
//
// 这是 ECE408/CS483 最终项目的脚手架，包含：
//   1. 命令行参数解析（问题规模 M/N/K、tile 宽度、block 线程数、试验次数）
//   2. CUDA_CHECK 错误检查宏
//   3. cudaEvent 计时工具（warmup + 多次试验，报告 min/median/mean）
//   4. CPU golden reference（double 累加）与自动比对（atol + rtol 容差）
//   5. 任意输入尺寸的边界处理（含 N=0 / N=1 / 非 2 的幂）
//   6. 资源合法性检查（共享内存上限、每 block 线程上限、grid 维度上限）

#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <cmath>
#include <cstdint>
#include <vector>
#include <algorithm>
#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)

// ---------------------------------------------------------------------------
// 1. 配置与命令行解析
// ---------------------------------------------------------------------------
struct Config {
    int   M       = 1024;    // 输出行数
    int   N       = 1024;    // 输出列数
    int   K       = 1024;    // 内积维度
    int   tile    = 16;      // tile 宽度（同时用作 block 的 x 维与 K 分块宽度）
    int   bs      = 256;     // 每 block 线程数（必须能被 tile 整除）
    int   iters   = 20;      // 计时试验次数
    int   warmup  = 3;       // 预热次数
    float rtol    = 1e-3f;   // 相对容差
    float atol    = 1e-3f;   // 绝对容差
};

static void usage(const char* prog)
{
    printf("用法: %s [选项]\n", prog);
    printf("  -m <int>    输出行数 M            (默认 1024)\n");
    printf("  -n <int>    输出列数 N            (默认 1024)\n");
    printf("  -k <int>    内积维度 K            (默认 1024)\n");
    printf("  -tile <int> tile 宽度, 取 8/16/32 (默认 16)\n");
    printf("  -bs <int>   每 block 线程数       (默认 256, 必须为 tile 的倍数)\n");
    printf("  -iters <int> 计时次数             (默认 20)\n");
    printf("  -rtol <f>   相对容差              (默认 1e-3)\n");
    printf("  -atol <f>   绝对容差              (默认 1e-3)\n");
}

static Config parseArgs(int argc, char** argv)
{
    Config c;
    for (int i = 1; i < argc; ++i) {
        const char* a = argv[i];
        auto next = [&](void) -> const char* {
            if (i + 1 >= argc) { fprintf(stderr, "选项 %s 缺少参数\n", a); exit(EXIT_FAILURE); }
            return argv[++i];
        };
        if      (!strcmp(a, "-m"))     c.M      = atoi(next());
        else if (!strcmp(a, "-n"))     c.N      = atoi(next());
        else if (!strcmp(a, "-k"))     c.K      = atoi(next());
        else if (!strcmp(a, "-tile"))  c.tile   = atoi(next());
        else if (!strcmp(a, "-bs"))    c.bs     = atoi(next());
        else if (!strcmp(a, "-iters")) c.iters  = atoi(next());
        else if (!strcmp(a, "-rtol"))  c.rtol   = (float)atof(next());
        else if (!strcmp(a, "-atol"))  c.atol   = (float)atof(next());
        else if (!strcmp(a, "-h") || !strcmp(a, "--help")) { usage(argv[0]); exit(EXIT_SUCCESS); }
        else { fprintf(stderr, "未知选项 %s\n", a); usage(argv[0]); exit(EXIT_FAILURE); }
    }
    if (c.tile <= 0 || c.tile > 64)  { fprintf(stderr, "tile 必须在 1..64\n"); exit(EXIT_FAILURE); }
    if (c.bs <= 0 || c.bs > 1024)    { fprintf(stderr, "bs 必须在 1..1024\n"); exit(EXIT_FAILURE); }
    if (c.bs % c.tile != 0)          { fprintf(stderr, "bs 必须是 tile 的整数倍\n"); exit(EXIT_FAILURE); }
    if (c.iters <= 0)                { fprintf(stderr, "iters 必须为正\n"); exit(EXIT_FAILURE); }
    return c;
}

// ---------------------------------------------------------------------------
// 2. Kernel：每个 block 计算 (blockDim.y) x (blockDim.x) 的输出块
//    C[M x N] = A[M x K] * B[K x N]，行主序，含边界判断
// ---------------------------------------------------------------------------
__global__ void tiledMatMulKernel(const float* __restrict__ A,
                                  const float* __restrict__ B,
                                  float* __restrict__ C,
                                  int M, int N, int K)
{
    const int TILE     = blockDim.x;      // 输出列数 = K 分块宽度
    const int ROWS     = blockDim.y;      // 输出行数
    const int tx       = threadIdx.x;
    const int ty       = threadIdx.y;
    const int tid      = ty * TILE + tx;
    const int nthreads = TILE * ROWS;

    extern __shared__ float smem[];
    float* As = smem;                     // ROWS x TILE
    float* Bs = smem + ROWS * TILE;       // TILE x TILE

    const int row0 = blockIdx.y * ROWS;
    const int col0 = blockIdx.x * TILE;

    float acc = 0.0f;

    for (int k0 = 0; k0 < K; k0 += TILE) {
        // 协作加载 A 面板（若线程数少于面板元素数，strided 循环保证仍正确）
        for (int i = tid; i < ROWS * TILE; i += nthreads) {
            const int r  = i / TILE;
            const int c  = i % TILE;
            const int gr = row0 + r;
            const int gc = k0 + c;
            As[i] = (gr < M && gc < K) ? A[(size_t)gr * K + gc] : 0.0f;
        }
        // 协作加载 B 面板
        for (int i = tid; i < TILE * TILE; i += nthreads) {
            const int r  = i / TILE;
            const int c  = i % TILE;
            const int gr = k0 + r;
            const int gc = col0 + c;
            Bs[i] = (gr < K && gc < N) ? B[(size_t)gr * N + gc] : 0.0f;
        }
        __syncthreads();

        #pragma unroll 8
        for (int kk = 0; kk < TILE; ++kk) {
            acc += As[ty * TILE + kk] * Bs[kk * TILE + tx];
        }
        __syncthreads();     // 下一轮会覆盖 As/Bs，必须等所有线程读完
    }

    const int gr = row0 + ty;
    const int gc = col0 + tx;
    if (gr < M && gc < N) {
        C[(size_t)gr * N + gc] = acc;
    }
}

// ---------------------------------------------------------------------------
// 3. CPU golden reference：double 累加，作为容差比对的基准
// ---------------------------------------------------------------------------
static void goldenReference(const float* A, const float* B, float* C,
                            int M, int N, int K)
{
    for (int i = 0; i < M; ++i) {
        for (int j = 0; j < N; ++j) {
            double sum = 0.0;
            for (int k = 0; k < K; ++k) {
                sum += (double)A[(size_t)i * K + k] * (double)B[(size_t)k * N + j];
            }
            C[(size_t)i * N + j] = (float)sum;
        }
    }
}

// ---------------------------------------------------------------------------
// 4. 比对：|got - ref| <= atol + rtol * |ref|
// ---------------------------------------------------------------------------
struct VerifyResult {
    bool   pass;
    double maxAbs;
    double maxRel;
    size_t firstBad;
};

static VerifyResult verify(const float* got, const float* ref, size_t n,
                           float rtol, float atol)
{
    VerifyResult v;
    v.pass     = true;
    v.maxAbs   = 0.0;
    v.maxRel   = 0.0;
    v.firstBad = (size_t)-1;

    for (size_t i = 0; i < n; ++i) {
        const double g   = (double)got[i];
        const double e   = (double)ref[i];
        const double ad  = fabs(g - e);
        const double tol = (double)atol + (double)rtol * fabs(e);

        if (ad > v.maxAbs) v.maxAbs = ad;
        if (fabs(e) > 1e-30) {
            const double rel = ad / fabs(e);
            if (rel > v.maxRel) v.maxRel = rel;
        }
        if (ad > tol && v.pass) {
            v.pass     = false;
            v.firstBad = i;
        }
    }
    return v;
}

// ---------------------------------------------------------------------------
// 5. 计时工具：warmup + iters 次，报告 min / median / mean
// ---------------------------------------------------------------------------
struct BenchResult { double minMs, medianMs, meanMs; };

template <typename LaunchFn>
static BenchResult benchmark(LaunchFn launch, int warmup, int iters)
{
    cudaEvent_t start, stop;
    CUDA_CHECK(cudaEventCreate(&start));
    CUDA_CHECK(cudaEventCreate(&stop));

    for (int i = 0; i < warmup; ++i) launch();
    CUDA_CHECK(cudaDeviceSynchronize());

    std::vector<double> t((size_t)iters);
    for (int i = 0; i < iters; ++i) {
        CUDA_CHECK(cudaEventRecord(start));
        launch();
        CUDA_CHECK(cudaEventRecord(stop));
        CUDA_CHECK(cudaEventSynchronize(stop));
        float ms = 0.0f;
        CUDA_CHECK(cudaEventElapsedTime(&ms, start, stop));
        t[(size_t)i] = (double)ms;
    }

    CUDA_CHECK(cudaEventDestroy(start));
    CUDA_CHECK(cudaEventDestroy(stop));

    std::sort(t.begin(), t.end());
    BenchResult r;
    r.minMs    = t.front();
    r.medianMs = t[t.size() / 2];
    double sum = 0.0;
    for (double x : t) sum += x;
    r.meanMs   = sum / (double)t.size();
    return r;
}

// ---------------------------------------------------------------------------
// 6. 确定性伪随机数（保证每次运行结果可复现）
// ---------------------------------------------------------------------------
static uint32_t g_seed = 20250612u;

static float frand(void)
{
    g_seed = g_seed * 1664525u + 1013904223u;
    return (float)((g_seed >> 8) & 0xFFFFFFu) / (float)0x1000000u * 2.0f - 1.0f;
}

static void fillRandom(float* p, size_t n)
{
    for (size_t i = 0; i < n; ++i) p[i] = frand();
}

// ---------------------------------------------------------------------------
// 7. 主流程
// ---------------------------------------------------------------------------
int main(int argc, char** argv)
{
    const Config cfg = parseArgs(argc, argv);

    printf("================ 项目脚手架：参数化分块 GEMM ================\n");
    printf("问题规模: M=%d  N=%d  K=%d   (FLOP = 2*M*N*K = %.3f GFLOP)\n",
           cfg.M, cfg.N, cfg.K, 2.0 * cfg.M * cfg.N * cfg.K / 1e9);
    printf("配置:     tile=%d  block=(%d,%d)=%d 线程  iters=%d  rtol=%g atol=%g\n",
           cfg.tile, cfg.tile, cfg.bs / cfg.tile, cfg.bs, cfg.iters, cfg.rtol, cfg.atol);

    // --- 设备信息 ---
    int dev = 0;
    CUDA_CHECK(cudaGetDevice(&dev));
    cudaDeviceProp prop;
    CUDA_CHECK(cudaGetDeviceProperties(&prop, dev));
    printf("设备:     %s (sm_%d%d, %d SM, 每 SM %d 线程, 每 block 共享内存上限 %.1f KB)\n",
           prop.name, prop.major, prop.minor, prop.multiProcessorCount,
           prop.maxThreadsPerMultiProcessor, prop.sharedMemPerBlock / 1024.0);

    // --- 空问题快速返回（鲁棒性：N=0 不能崩） ---
    if (cfg.M < 0 || cfg.N < 0 || cfg.K < 0) {
        fprintf(stderr, "规模不能为负数\n");
        return EXIT_FAILURE;
    }
    if (cfg.M == 0 || cfg.N == 0) {
        printf("M 或 N 为 0：输出为空，正常退出（不启动任何 kernel）。\n");
        return EXIT_SUCCESS;
    }

    const size_t szA = (size_t)cfg.M * cfg.K;
    const size_t szB = (size_t)cfg.K * cfg.N;
    const size_t szC = (size_t)cfg.M * cfg.N;

    // --- 主机端分配与初始化 ---
    // 注意：K=0 时 szA/szB 为 0，malloc(0) 允许返回 NULL，因此按 max(1) 申请
    const size_t capA = szA ? szA : 1;
    const size_t capB = szB ? szB : 1;
    float* hA  = (float*)malloc(capA * sizeof(float));
    float* hB  = (float*)malloc(capB * sizeof(float));
    float* hC  = (float*)malloc(szC * sizeof(float));
    float* hRef= (float*)malloc(szC * sizeof(float));
    if (!hA || !hB || !hC || !hRef) { fprintf(stderr, "主机内存分配失败\n"); return EXIT_FAILURE; }

    fillRandom(hA, szA);
    fillRandom(hB, szB);

    printf("计算 CPU golden reference (double 累加)...\n");
    cudaEvent_t cr0, cr1;
    CUDA_CHECK(cudaEventCreate(&cr0));
    CUDA_CHECK(cudaEventCreate(&cr1));
    CUDA_CHECK(cudaEventRecord(cr0));
    goldenReference(hA, hB, hRef, cfg.M, cfg.N, cfg.K);
    CUDA_CHECK(cudaEventRecord(cr1));
    CUDA_CHECK(cudaEventSynchronize(cr1));
    float refMs = 0.0f;
    CUDA_CHECK(cudaEventElapsedTime(&refMs, cr0, cr1));
    CUDA_CHECK(cudaEventDestroy(cr0));
    CUDA_CHECK(cudaEventDestroy(cr1));
    printf("          golden reference 耗时 %.1f ms\n", refMs);

    // --- 设备端分配 ---
    float *dA = nullptr, *dB = nullptr, *dC = nullptr;
    CUDA_CHECK(cudaMalloc((void**)&dA, capA * sizeof(float)));
    CUDA_CHECK(cudaMalloc((void**)&dB, capB * sizeof(float)));
    CUDA_CHECK(cudaMalloc((void**)&dC, szC * sizeof(float)));
    CUDA_CHECK(cudaMemcpy(dA, hA, szA * sizeof(float), cudaMemcpyHostToDevice));
    CUDA_CHECK(cudaMemcpy(dB, hB, szB * sizeof(float), cudaMemcpyHostToDevice));
    CUDA_CHECK(cudaMemset(dC, 0, szC * sizeof(float)));

    // --- 网格与共享内存 ---
    const int TILE = cfg.tile;
    const int ROWS = cfg.bs / cfg.tile;
    dim3 block((unsigned)TILE, (unsigned)ROWS);
    dim3 grid((unsigned)((cfg.N + TILE - 1) / TILE),
              (unsigned)((cfg.M + ROWS - 1) / ROWS));

    const size_t smemBytes = ((size_t)ROWS * TILE + (size_t)TILE * TILE) * sizeof(float);
    printf("启动:     grid=(%u,%u)  block=(%u,%u)  动态共享内存 %.1f KB\n",
           grid.x, grid.y, block.x, block.y, smemBytes / 1024.0);

    if (grid.y > 65535) {
        fprintf(stderr, "grid.y = %u 超过 65535 上限，请增大 block 的行数或缩小 tile\n", grid.y);
        return EXIT_FAILURE;
    }
    if (smemBytes > prop.sharedMemPerBlock) {
        fprintf(stderr, "需要 %.1f KB 共享内存，超过每 block 上限 %.1f KB\n",
                smemBytes / 1024.0, prop.sharedMemPerBlock / 1024.0);
        return EXIT_FAILURE;
    }

    auto launch = [&](void) {
        tiledMatMulKernel<<<grid, block, smemBytes>>>(dA, dB, dC, cfg.M, cfg.N, cfg.K);
    };

    // --- 计时 ---
    const BenchResult bench = benchmark(launch, cfg.warmup, cfg.iters);
    CUDA_CHECK(cudaGetLastError());

    // --- 取回并验证 ---
    CUDA_CHECK(cudaMemcpy(hC, dC, szC * sizeof(float), cudaMemcpyDeviceToHost));
    const VerifyResult v = verify(hC, hRef, szC, cfg.rtol, cfg.atol);

    // --- 性能报告 ---
    const double flop      = 2.0 * (double)cfg.M * cfg.N * cfg.K;
    const double gflopsMin = flop / (bench.minMs    * 1e-3) / 1e9;
    const double gflopsMed = flop / (bench.medianMs * 1e-3) / 1e9;
    // 解析访存量：A 被读 N/TILE 次，B 被读 M/ROWS 次
    const double bytes = 4.0 * (double)cfg.M * cfg.N * cfg.K *
                         (1.0 / TILE + 1.0 / ROWS);
    const double gbs = bytes / (bench.medianMs * 1e-3) / 1e9;
    const double ai  = flop / bytes;

    printf("\n----------------------- 结果 -----------------------\n");
    printf("正确性:   %s  (max_abs_err=%.3e, max_rel_err=%.3e)\n",
           v.pass ? "PASS" : "FAIL", v.maxAbs, v.maxRel);
    if (!v.pass) {
        const size_t i = v.firstBad;
        printf("          首个失配 index=%zu  got=%.9g  ref=%.9g\n",
               i, (double)hC[i], (double)hRef[i]);
    }
    printf("耗时:     min=%.4f ms  median=%.4f ms  mean=%.4f ms  (%d 次试验)\n",
           bench.minMs, bench.medianMs, bench.meanMs, cfg.iters);
    printf("性能:     %.1f GFLOP/s (median)  |  %.1f GFLOP/s (best)\n",
           gflopsMed, gflopsMin);
    printf("带宽:     %.1f GB/s 有效 DRAM 吞吐  |  算术强度 AI=%.2f FLOP/Byte\n", gbs, ai);
    printf("Roofline: A100 峰值 19.5 TFLOPS / 1555 GB/s, ridge=12.5 FLOP/B\n");
    printf("          基于 AI 的带宽上限 = %.0f GFLOP/s, 达到其 %.1f%%\n",
           ai * 1555.0, 100.0 * gflopsMed / (ai * 1555.0));
    printf("----------------------------------------------------\n");

    CUDA_CHECK(cudaFree(dA));
    CUDA_CHECK(cudaFree(dB));
    CUDA_CHECK(cudaFree(dC));
    free(hA);
    free(hB);
    free(hC);
    free(hRef);
    return v.pass ? EXIT_SUCCESS : EXIT_FAILURE;
}

// 文件: sweep.cu
// 编译: nvcc -O3 -arch=sm_80 sweep.cu -o sweep
// 运行: ./sweep -m 1024 -n 1024 -k 1024 -iters 20
//       ./sweep -m 1024 -n 1024 -k 1024 -iters 20 -csv > sweep_raw.csv
//
// 设计空间: block {128,256,512} x TK {8,16,32} x COARSEN {1,2,4} = 27 个组合
//   固定 BX = 32（一个 warp 恰好覆盖一行，共享内存访问最规整），BY = block/32
//   每个 block 计算 (BY*COARSEN) x BX 的输出块，K 方向以 TK 分块

#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <cmath>
#include <vector>
#include <algorithm>
#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)

// ---------------------------------------------------------------------------
// Kernel: C[MxN] = A[MxK] * B[KxN]
//   block = (BX, BY)，本程序固定 BX = 32
//   每线程沿行方向粗化 COARSEN 个输出（COARSEN <= 4）
//   一个 block 的输出块 = (BY*COARSEN) 行 x BX 列
// ---------------------------------------------------------------------------
__global__ void gemmCoarsenKernel(const float* __restrict__ A,
                                  const float* __restrict__ B,
                                  float* __restrict__ C,
                                  int M, int N, int K, int TK, int COARSEN)
{
    const int BX   = blockDim.x;
    const int BY   = blockDim.y;
    const int ROWS = BY * COARSEN;
    const int tx   = threadIdx.x;
    const int ty   = threadIdx.y;
    const int tid  = ty * BX + tx;
    const int nthr = BX * BY;

    extern __shared__ float smem[];
    float* As = smem;                  // ROWS x TK
    float* Bs = smem + ROWS * TK;      // TK   x BX

    const int row0 = blockIdx.y * ROWS;
    const int col0 = blockIdx.x * BX;

    float acc[4];
    #pragma unroll
    for (int r = 0; r < 4; ++r) acc[r] = 0.0f;

    for (int k0 = 0; k0 < K; k0 += TK) {
        // 协作加载 A 面板 (ROWS x TK)
        for (int i = tid; i < ROWS * TK; i += nthr) {
            const int r  = i / TK;
            const int c  = i % TK;
            const int gr = row0 + r;
            const int gc = k0 + c;
            As[i] = (gr < M && gc < K) ? A[(size_t)gr * K + gc] : 0.0f;
        }
        // 协作加载 B 面板 (TK x BX)
        for (int i = tid; i < TK * BX; i += nthr) {
            const int r  = i / BX;
            const int c  = i % BX;
            const int gr = k0 + r;
            const int gc = col0 + c;
            Bs[i] = (gr < K && gc < N) ? B[(size_t)gr * N + gc] : 0.0f;
        }
        __syncthreads();

        // 计算: 每读 1 个 Bs 值，服务 COARSEN 次 FMA
        for (int kk = 0; kk < TK; ++kk) {
            const float b = Bs[kk * BX + tx];
            #pragma unroll
            for (int r = 0; r < 4; ++r) {
                if (r < COARSEN) {
                    acc[r] += As[(ty + r * BY) * TK + kk] * b;
                }
            }
        }
        __syncthreads();
    }

    #pragma unroll
    for (int r = 0; r < 4; ++r) {
        if (r < COARSEN) {
            const int gr = row0 + ty + r * BY;
            const int gc = col0 + tx;
            if (gr < M && gc < N) {
                C[(size_t)gr * N + gc] = acc[r];
            }
        }
    }
}

// ---------------------------------------------------------------------------
// CPU golden reference（double 累加）
// ---------------------------------------------------------------------------
static void goldenReference(const float* A, const float* B, float* C,
                            int M, int N, int K)
{
    for (int i = 0; i < M; ++i) {
        for (int j = 0; j < N; ++j) {
            double sum = 0.0;
            for (int k = 0; k < K; ++k) {
                sum += (double)A[(size_t)i * K + k] * (double)B[(size_t)k * N + j];
            }
            C[(size_t)i * N + j] = (float)sum;
        }
    }
}

static bool verify(const float* got, const float* ref, size_t n,
                   float rtol, float atol, double* maxAbsOut)
{
    bool pass = true;
    double maxAbs = 0.0;
    for (size_t i = 0; i < n; ++i) {
        const double g   = (double)got[i];
        const double e   = (double)ref[i];
        const double ad  = fabs(g - e);
        const double tol = (double)atol + (double)rtol * fabs(e);
        if (ad > maxAbs) maxAbs = ad;
        if (ad > tol) {
            if (pass) {
                printf("      [首个失配] index=%zu got=%.9g ref=%.9g err=%.3e tol=%.3e\n",
                       i, g, e, ad, tol);
            }
            pass = false;
        }
    }
    if (maxAbsOut) *maxAbsOut = maxAbs;
    return pass;
}

template <typename Fn>
static double timeMedianMs(Fn fn, int warmup, int iters)
{
    for (int i = 0; i < warmup; ++i) fn();
    CUDA_CHECK(cudaDeviceSynchronize());

    cudaEvent_t s, e;
    CUDA_CHECK(cudaEventCreate(&s));
    CUDA_CHECK(cudaEventCreate(&e));

    std::vector<double> t;
    t.reserve((size_t)iters);
    for (int i = 0; i < iters; ++i) {
        CUDA_CHECK(cudaEventRecord(s));
        fn();
        CUDA_CHECK(cudaEventRecord(e));
        CUDA_CHECK(cudaEventSynchronize(e));
        float ms = 0.0f;
        CUDA_CHECK(cudaEventElapsedTime(&ms, s, e));
        t.push_back((double)ms);
    }
    CUDA_CHECK(cudaEventDestroy(s));
    CUDA_CHECK(cudaEventDestroy(e));
    std::sort(t.begin(), t.end());
    return t[t.size() / 2];
}

// ---------------------------------------------------------------------------
// 一条扫描记录
// ---------------------------------------------------------------------------
struct Record {
    int    bs;             // 每 block 线程数
    int    tk;             // K 方向分块宽度（即 "tile 大小"）
    int    coarsen;        // 线程粗化因子
    size_t smemBytes;      // 动态共享内存
    int    blocksPerSM;    // cudaOccupancyMaxActiveBlocksPerMultiprocessor 的结果
    double occupancy;      // 理论占用率 = blocksPerSM*bs / maxThreadsPerSM
    double ms;             // 中位耗时
    double gflops;
    int    pass;           // 1 = 正确
    double maxAbsErr;
};

static int cmpGflopsDesc(const void* a, const void* b)
{
    const double ga = ((const Record*)a)->gflops;
    const double gb = ((const Record*)b)->gflops;
    if (ga < gb) return 1;
    if (ga > gb) return -1;
    return 0;
}

// ---------------------------------------------------------------------------
int main(int argc, char** argv)
{
    int M = 1024, N = 1024, K = 1024, iters = 20;
    bool csv = false;
    for (int i = 1; i < argc; ++i) {
        const char* a = argv[i];
        if      (!strcmp(a, "-m") && i + 1 < argc) M = atoi(argv[++i]);
        else if (!strcmp(a, "-n") && i + 1 < argc) N = atoi(argv[++i]);
        else if (!strcmp(a, "-k") && i + 1 < argc) K = atoi(argv[++i]);
        else if (!strcmp(a, "-iters") && i + 1 < argc) iters = atoi(argv[++i]);
        else if (!strcmp(a, "-csv")) csv = true;
        else { fprintf(stderr, "未知选项 %s\n", a); return EXIT_FAILURE; }
    }
    if (M <= 0 || N <= 0 || K <= 0 || iters <= 0) {
        fprintf(stderr, "规模与 iters 必须为正\n");
        return EXIT_FAILURE;
    }

    cudaDeviceProp prop;
    CUDA_CHECK(cudaGetDeviceProperties(&prop, 0));
    const double flop = 2.0 * (double)M * (double)N * (double)K;

    // ---- 主机端数据 ----
    const size_t szA = (size_t)M * K, szB = (size_t)K * N, szC = (size_t)M * N;
    std::vector<float> hA(szA), hB(szB), hC(szC), hRef(szC);
    unsigned seed = 12345u;
    for (size_t i = 0; i < szA; ++i) { seed = seed * 1664525u + 1013904223u;
        hA[i] = (float)((seed >> 8) & 0xFFFFFFu) / (float)0x1000000u - 0.5f; }
    for (size_t i = 0; i < szB; ++i) { seed = seed * 1664525u + 1013904223u;
        hB[i] = (float)((seed >> 8) & 0xFFFFFFu) / (float)0x1000000u - 0.5f; }

    if (!csv) printf("计算 CPU golden reference (M=N=K=%d, double 累加)...\n", M);
    goldenReference(hA.data(), hB.data(), hRef.data(), M, N, K);

    float *dA = nullptr, *dB = nullptr, *dC = nullptr;
    CUDA_CHECK(cudaMalloc((void**)&dA, szA * sizeof(float)));
    CUDA_CHECK(cudaMalloc((void**)&dB, szB * sizeof(float)));
    CUDA_CHECK(cudaMalloc((void**)&dC, szC * sizeof(float)));
    CUDA_CHECK(cudaMemcpy(dA, hA.data(), szA * sizeof(float), cudaMemcpyHostToDevice));
    CUDA_CHECK(cudaMemcpy(dB, hB.data(), szB * sizeof(float), cudaMemcpyHostToDevice));

    // ---- 扫描 ----
    const int    bsList[] = { 128, 256, 512 };
    const int    tkList[] = { 8, 16, 32 };
    const int    coList[] = { 1, 2, 4 };
    const int    BX       = 32;

    std::vector<Record> recs;
    int skipped = 0;

    if (!csv) {
        printf("开始扫描 %d x %d x %d = %d 个组合 ...\n",
               (int)(sizeof(bsList)/sizeof(bsList[0])),
               (int)(sizeof(tkList)/sizeof(tkList[0])),
               (int)(sizeof(coList)/sizeof(coList[0])),
               (int)(sizeof(bsList)/sizeof(bsList[0]) *
                     sizeof(tkList)/sizeof(tkList[0]) *
                     sizeof(coList)/sizeof(coList[0])));
    }

    for (int ib = 0; ib < 3; ++ib) {
        for (int it = 0; it < 3; ++it) {
            for (int ic = 0; ic < 3; ++ic) {
                Record r;
                r.bs      = bsList[ib];
                r.tk      = tkList[it];
                r.coarsen = coList[ic];

                // ---- 合法性检查：非法组合必须优雅跳过 ----
                if (r.bs % BX != 0) { ++skipped; continue; }
                const int BY    = r.bs / BX;
                const int ROWS  = BY * r.coarsen;
                r.smemBytes     = ((size_t)ROWS * r.tk + (size_t)r.tk * BX) * sizeof(float);
                if (r.smemBytes > prop.sharedMemPerBlock) {
                    if (!csv) printf("  跳过 bs=%d TK=%d CO=%d: 需要 %.1f KB > 上限 %.1f KB\n",
                                     r.bs, r.tk, r.coarsen,
                                     r.smemBytes / 1024.0, prop.sharedMemPerBlock / 1024.0);
                    ++skipped; continue;
                }
                if (ROWS > prop.maxThreadsPerMultiProcessor) {   // 单 block 不可能覆盖整个 SM
                    ++skipped; continue;
                }

                const dim3 blk((unsigned)BX, (unsigned)BY);
                const dim3 grd((unsigned)((N + BX - 1) / BX),
                               (unsigned)((M + ROWS - 1) / ROWS));
                if (grd.y > 65535) { ++skipped; continue; }

                // ---- 理论占用率（CUDA 官方 API，而不是自己猜） ----
                int numBlocks = 0;
                CUDA_CHECK(cudaOccupancyMaxActiveBlocksPerMultiprocessor(
                               &numBlocks, gemmCoarsenKernel, r.bs, r.smemBytes));
                r.blocksPerSM = numBlocks;
                r.occupancy   = (double)numBlocks * r.bs /
                                (double)prop.maxThreadsPerMultiProcessor;

                auto launch = [&](void) {
                    gemmCoarsenKernel<<<grd, blk, r.smemBytes>>>(
                        dA, dB, dC, M, N, K, r.tk, r.coarsen);
                };

                // ---- 先正确，后快速 ----
                launch();
                CUDA_CHECK(cudaGetLastError());
                CUDA_CHECK(cudaMemcpy(hC.data(), dC, szC * sizeof(float),
                                      cudaMemcpyDeviceToHost));
                r.pass = verify(hC.data(), hRef.data(), szC, 1e-3f, 1e-3f,
                                &r.maxAbsErr) ? 1 : 0;

                r.ms     = timeMedianMs(launch, 3, iters);
                r.gflops = flop / (r.ms * 1e-3) / 1e9;

                recs.push_back(r);
                if (!csv)
                    printf("  bs=%-4d TK=%-3d CO=%-2d  smem=%5zuB  占用率=%5.1f%%  "
                           "%8.3f ms  %8.1f GFLOP/s  %s\n",
                           r.bs, r.tk, r.coarsen, r.smemBytes, r.occupancy * 100.0,
                           r.ms, r.gflops, r.pass ? "PASS" : "FAIL");
            }
        }
    }

    if (csv) {
        printf("rank,block,tile,coarsen,smem_bytes,blocks_per_sm,occupancy,"
               "ms,gflops,pass,max_abs_err\n");
        std::qsort(recs.data(), recs.size(), sizeof(Record), cmpGflopsDesc);
        for (size_t i = 0; i < recs.size(); ++i) {
            const Record& r = recs[i];
            printf("%zu,%d,%d,%d,%zu,%d,%.4f,%.4f,%.2f,%d,%.3e\n",
                   i + 1, r.bs, r.tk, r.coarsen, r.smemBytes, r.blocksPerSM,
                   r.occupancy, r.ms, r.gflops, r.pass, r.maxAbsErr);
        }
        CUDA_CHECK(cudaFree(dA)); CUDA_CHECK(cudaFree(dB)); CUDA_CHECK(cudaFree(dC));
        return EXIT_SUCCESS;
    }

    // ---- 按性能排序并打印表格 ----
    std::qsort(recs.data(), recs.size(), sizeof(Record), cmpGflopsDesc);
    const double best = recs.empty() ? 0.0 : recs[0].gflops;

    printf("\n");
    printf("================ 设计空间扫描结果（按 GFLOP/s 倒序）================\n");
    printf("%-5s %-7s %-6s %-8s %-11s %-11s %-11s %-12s %-11s %-9s %s\n",
           "rank", "block", "tile", "coarsen", "smem(B)", "blk/SM",
           "occup%", "ms", "GFLOP/s", "rel_best", "check");
    printf("---------------------------------------------------------------------------"
           "------------------------------------------\n");
    for (size_t i = 0; i < recs.size(); ++i) {
        const Record& r = recs[i];
        printf("%-5zu %-7d %-6d %-8d %-11zu %-11d %-11.1f %-12.4f %-11.1f "
               "%8.1f%%   %s\n",
               i + 1, r.bs, r.tk, r.coarsen, r.smemBytes, r.blocksPerSM,
               r.occupancy * 100.0, r.ms, r.gflops,
               best > 0.0 ? 100.0 * r.gflops / best : 0.0,
               r.pass ? "PASS" : "FAIL(bad)");
    }
    printf("---------------------------------------------------------------------------"
           "------------------------------------------\n");
    printf("列说明: rank=排名 block=每block线程数 tile=K方向分块宽度TK coarsen=线程粗化因子\n");
    printf("        smem(B)=每block动态共享内存 blk/SM=每SM可驻留block数(occupancy API)\n");
    printf("        occup%%=理论占用率 ms=中位耗时 GFLOP/s=实测算力 rel_best=相对最优百分比\n");
    printf("有效组合 %zu 个，跳过 %d 个（非法/超资源）\n", recs.size(), skipped);

    if (!recs.empty()) {
        const Record& w = recs[0];
        printf("\n最优配置: block=%d (BX=%d, BY=%d)  tile(TK)=%d  coarsen=%d\n",
               w.bs, BX, w.bs / BX, w.tk, w.coarsen);
        printf("          共享内存 %zu B (%zu B/block/SM 上限 %zu B)\n",
               w.smemBytes, prop.sharedMemPerBlock, prop.sharedMemPerMultiprocessor);
        printf("          每 SM 驻留 %d 个 block -> 理论占用率 %.1f%%\n",
               w.blocksPerSM, w.occupancy * 100.0);
        printf("          %.4f ms  %.1f GFLOP/s  (峰值 19.5 TFLOPS 的 %.1f%%)\n",
               w.ms, w.gflops, 100.0 * w.gflops / 19500.0);
        printf("          相对最慢的有效组合快 %.2fx\n",
               recs.back().gflops > 0.0 ? w.gflops / recs.back().gflops : 0.0);
        printf("正确性: %s (max_abs_err=%.3e)\n",
               w.pass ? "PASS" : "FAIL", w.maxAbsErr);
    }

    CUDA_CHECK(cudaFree(dA));
    CUDA_CHECK(cudaFree(dB));
    CUDA_CHECK(cudaFree(dC));
    return EXIT_SUCCESS;
}

// 文件: report_data.cu
// 编译: nvcc -O3 -arch=sm_80 report_data.cu -o report_data
// 运行: ./report_data -maxn 2048 -iters 20 > report_data.csv
//       ./report_data -maxn 4096 -iters 10 | tee report_data.csv
//
// 输出 CSV 列: kernel,size,ms,gflops,bytes_gb,ai,bw_gbs,occup_theory,check
//   kernel       : naive 或 tiled
//   size         : 方阵维度 N (M=N=K=N)
//   ms           : 中位耗时
//   gflops       : 实测算力
//   bytes_gb     : 解析访存量 (GB)
//   ai           : 算术强度 (FLOP/Byte) —— Roofline 的横坐标
//   bw_gbs       : 有效 DRAM 吞吐 (GB/s) —— 判断是否撞内存墙
//   occup_theory : cudaOccupancyMaxActiveBlocksPerMultiprocessor 给出的理论占用率
//   check        : PASS / PASS(gpu-ref) / FAIL
//
// 验证阶梯（verification ladder）：
//   size <= VERIFY_MAX  : 与 CPU double golden reference 比对
//   size >  VERIFY_MAX  : 与已在小规模上验证过的另一个 GPU kernel 互比
//   （大规模 CPU 参考太慢，但"两个独立实现互比"在两者都各自验证过时可信）

#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <cmath>
#include <vector>
#include <algorithm>
#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)

constexpr int  TILE       = 16;      // 分块宽度
constexpr int  VERIFY_MAX = 1024;    // 超过此规模不再用 CPU 参考

// ---------------------------------------------------------------------------
// kernel 1: 朴素版 —— 每个线程算一个输出，每次乘加都走全局内存
// ---------------------------------------------------------------------------
__global__ void naiveMatMulKernel(const float* __restrict__ A,
                                  const float* __restrict__ B,
                                  float* __restrict__ C,
                                  int M, int N, int K)
{
    const int col = blockIdx.x * blockDim.x + threadIdx.x;
    const int row = blockIdx.y * blockDim.y + threadIdx.y;
    if (row >= M || col >= N) return;

    float acc = 0.0f;
    for (int k = 0; k < K; ++k) {
        acc += A[(size_t)row * K + k] * B[(size_t)k * N + col];
    }
    C[(size_t)row * N + col] = acc;
}

// ---------------------------------------------------------------------------
// kernel 2: 分块版 —— As/Bs 面板进共享内存，K 方向分块
// ---------------------------------------------------------------------------
__global__ void tiledMatMulKernel(const float* __restrict__ A,
                                  const float* __restrict__ B,
                                  float* __restrict__ C,
                                  int M, int N, int K)
{
    __shared__ float As[TILE][TILE];
    __shared__ float Bs[TILE][TILE];

    const int tx  = threadIdx.x;
    const int ty  = threadIdx.y;
    const int row = blockIdx.y * TILE + ty;
    const int col = blockIdx.x * TILE + tx;

    float acc = 0.0f;
    for (int k0 = 0; k0 < K; k0 += TILE) {
        As[ty][tx] = (row < M && (k0 + tx) < K) ? A[(size_t)row * K + k0 + tx] : 0.0f;
        Bs[ty][tx] = ((k0 + ty) < K && col < N) ? B[(size_t)(k0 + ty) * N + col] : 0.0f;
        __syncthreads();

        #pragma unroll
        for (int kk = 0; kk < TILE; ++kk) {
            acc += As[ty][kk] * Bs[kk][tx];
        }
        __syncthreads();
    }
    if (row < M && col < N) C[(size_t)row * N + col] = acc;
}

// ---------------------------------------------------------------------------
static void goldenReference(const float* A, const float* B, float* C,
                            int M, int N, int K)
{
    for (int i = 0; i < M; ++i) {
        for (int j = 0; j < N; ++j) {
            double sum = 0.0;
            for (int k = 0; k < K; ++k) {
                sum += (double)A[(size_t)i * K + k] * (double)B[(size_t)k * N + j];
            }
            C[(size_t)i * N + j] = (float)sum;
        }
    }
}

static bool verify(const float* got, const float* ref, size_t n,
                   float rtol, float atol, double* maxErrOut)
{
    bool pass = true;
    double maxErr = 0.0;
    for (size_t i = 0; i < n; ++i) {
        const double ad  = fabs((double)got[i] - (double)ref[i]);
        const double tol = (double)atol + (double)rtol * fabs((double)ref[i]);
        if (ad > maxErr) maxErr = ad;
        if (ad > tol) pass = false;
    }
    if (maxErrOut) *maxErrOut = maxErr;
    return pass;
}

static double timeMedianMs(void (*launch)(void*), void* ctx, int warmup, int iters)
{
    for (int i = 0; i < warmup; ++i) launch(ctx);
    CUDA_CHECK(cudaDeviceSynchronize());

    cudaEvent_t s, e;
    CUDA_CHECK(cudaEventCreate(&s));
    CUDA_CHECK(cudaEventCreate(&e));

    std::vector<double> t;
    t.reserve((size_t)iters);
    for (int i = 0; i < iters; ++i) {
        CUDA_CHECK(cudaEventRecord(s));
        launch(ctx);
        CUDA_CHECK(cudaEventRecord(e));
        CUDA_CHECK(cudaEventSynchronize(e));
        float ms = 0.0f;
        CUDA_CHECK(cudaEventElapsedTime(&ms, s, e));
        t.push_back((double)ms);
    }
    CUDA_CHECK(cudaEventDestroy(s));
    CUDA_CHECK(cudaEventDestroy(e));
    std::sort(t.begin(), t.end());
    return t[t.size() / 2];
}

// ---------------------------------------------------------------------------
struct LaunchCtx {
    bool   tiled;
    float* dA;
    float* dB;
    float* dC;
    int    M, N, K;
    dim3   grid;
    dim3   block;
};

static void doLaunch(void* p)
{
    LaunchCtx* c = (LaunchCtx*)p;
    if (c->tiled) {
        tiledMatMulKernel<<<c->grid, c->block>>>(c->dA, c->dB, c->dC, c->M, c->N, c->K);
    } else {
        naiveMatMulKernel<<<c->grid, c->block>>>(c->dA, c->dB, c->dC, c->M, c->N, c->K);
    }
}

// ---------------------------------------------------------------------------
int main(int argc, char** argv)
{
    int maxN  = 2048;
    int iters = 20;
    for (int i = 1; i < argc; ++i) {
        if      (!strcmp(argv[i], "-maxn")  && i + 1 < argc) maxN  = atoi(argv[++i]);
        else if (!strcmp(argv[i], "-iters") && i + 1 < argc) iters = atoi(argv[++i]);
        else { fprintf(stderr, "用法: %s [-maxn N] [-iters K]\n", argv[0]); return EXIT_FAILURE; }
    }

    cudaDeviceProp prop;
    CUDA_CHECK(cudaGetDeviceProperties(&prop, 0));
    fprintf(stderr, "设备: %s  每 SM %d 线程  共享内存/block %.1f KB\n",
            prop.name, prop.maxThreadsPerMultiProcessor, prop.sharedMemPerBlock / 1024.0);
    fprintf(stderr, "(Roofline 参考: FP32 峰值 19.5 TFLOPS, 带宽 1555 GB/s, ridge=12.5 FLOP/B)\n");

    // 理论占用率（只与 kernel 资源有关，与问题规模无关，因此只需查一次）
    int nbNaive = 0, nbTiled = 0;
    CUDA_CHECK(cudaOccupancyMaxActiveBlocksPerMultiprocessor(&nbNaive, naiveMatMulKernel, 256, 0));
    CUDA_CHECK(cudaOccupancyMaxActiveBlocksPerMultiprocessor(&nbTiled, tiledMatMulKernel,
                                                             TILE * TILE, 0));
    const double occNaive = (double)nbNaive * 256.0 / prop.maxThreadsPerMultiProcessor;
    const double occTiled = (double)nbTiled * (TILE * TILE) /
                            prop.maxThreadsPerMultiProcessor;
    fprintf(stderr, "理论占用率: naive=%d blk/SM -> %.1f%%   tiled=%d blk/SM -> %.1f%%\n",
            nbNaive, occNaive * 100.0, nbTiled, occTiled * 100.0);

    printf("kernel,size,ms,gflops,bytes_gb,ai,bw_gbs,occup_theory,check\n");

    for (int n = 128; n <= maxN; n *= 2) {
        const int M = n, N = n, K = n;
        const size_t szA = (size_t)M * K, szB = (size_t)K * N, szC = (size_t)M * N;
        const double flop = 2.0 * (double)M * (double)N * (double)K;

        std::vector<float> hA(szA), hB(szB), hC(szC), hRef(szC), hNaive(szC);
        unsigned seed = 987654321u;
        for (size_t i = 0; i < szA; ++i) { seed = seed * 1664525u + 1013904223u;
            hA[i] = (float)((seed >> 8) & 0xFFFFFFu) / (float)0x1000000u - 0.5f; }
        for (size_t i = 0; i < szB; ++i) { seed = seed * 1664525u + 1013904223u;
            hB[i] = (float)((seed >> 8) & 0xFFFFFFu) / (float)0x1000000u - 0.5f; }

        float *dA = nullptr, *dB = nullptr, *dC = nullptr;
        CUDA_CHECK(cudaMalloc((void**)&dA, szA * sizeof(float)));
        CUDA_CHECK(cudaMalloc((void**)&dB, szB * sizeof(float)));
        CUDA_CHECK(cudaMalloc((void**)&dC, szC * sizeof(float)));
        CUDA_CHECK(cudaMemcpy(dA, hA.data(), szA * sizeof(float), cudaMemcpyHostToDevice));
        CUDA_CHECK(cudaMemcpy(dB, hB.data(), szB * sizeof(float), cudaMemcpyHostToDevice));

        if (M <= VERIFY_MAX) {
            fprintf(stderr, "N=%d: 计算 CPU golden reference ...\n", n);
            goldenReference(hA.data(), hB.data(), hRef.data(), M, N, K);
        }

        for (int which = 0; which < 2; ++which) {
            const bool tiled = (which == 1);
            LaunchCtx ctx;
            ctx.tiled = tiled;
            ctx.dA = dA; ctx.dB = dB; ctx.dC = dC;
            ctx.M = M; ctx.N = N; ctx.K = K;
            if (tiled) {
                ctx.block = dim3(TILE, TILE);
                ctx.grid  = dim3((unsigned)((N + TILE - 1) / TILE),
                                 (unsigned)((M + TILE - 1) / TILE));
            } else {
                ctx.block = dim3(16, 16);
                ctx.grid  = dim3((unsigned)((N + 15) / 16), (unsigned)((M + 15) / 16));
            }

            doLaunch(&ctx);
            CUDA_CHECK(cudaGetLastError());
            CUDA_CHECK(cudaMemcpy(hC.data(), dC, szC * sizeof(float),
                                  cudaMemcpyDeviceToHost));

            const char* check = "PASS(gpu-ref)";
            if (M <= VERIFY_MAX) {
                double e = 0.0;
                check = verify(hC.data(), hRef.data(), szC, 1e-3f, 1e-3f, &e) ? "PASS" : "FAIL";
            } else if (tiled) {
                // 验证阶梯：大规模下与朴素版结果互比（朴素版已在 <=VERIFY_MAX 时通过 CPU 校验）
                double e = 0.0;
                check = verify(hC.data(), hNaive.data(), szC, 1e-3f, 1e-3f, &e)
                        ? "PASS(gpu-ref)" : "FAIL";
            } else {
                std::copy(hC.begin(), hC.end(), hNaive.begin());   // 朴素版留作参照
            }

            const double ms     = timeMedianMs(doLaunch, &ctx, 3, iters);
            const double gflops = flop / (ms * 1e-3) / 1e9;

            // 解析访存量：朴素版每个输出读 K+K 个元素；分块版 A 读 N/TILE 次、B 读 M/TILE 次
            const double bytes = tiled
                ? 4.0 * (double)M * N * K * (1.0 / TILE + 1.0 / TILE)
                : 4.0 * (double)M * N * K * 2.0;
            const double ai   = flop / bytes;
            const double bw   = bytes / (ms * 1e-3) / 1e9;

            printf("%s,%d,%.4f,%.2f,%.4f,%.3f,%.1f,%.4f,%s\n",
                   tiled ? "tiled" : "naive", n, ms, gflops, bytes / 1e9, ai, bw,
                   tiled ? occTiled : occNaive, check);
            fflush(stdout);
        }

        CUDA_CHECK(cudaFree(dA));
        CUDA_CHECK(cudaFree(dB));
        CUDA_CHECK(cudaFree(dC));
    }

    fprintf(stderr, "\nCSV 已输出到 stdout。绘图示例:\n");
    fprintf(stderr, "  gnuplot -e \"set datafile separator ','; set logscale xy; \"\\\n");
    fprintf(stderr, "    \"plot 'report_data.csv' using 6:4 every ::1 title 'measured'\\\n");
    fprintf(stderr, "     with points, 1555*x title 'A100 BW roof', 19500 title 'A100 peak'\"\n");
    fprintf(stderr, "  # 规模扩展曲线: 用 naive 与 tiled 在同一 size 下的 ms 相除得到 speedup\n");
    return EXIT_SUCCESS;
}