Lecture 9: 并行模式七 —— 稀疏矩阵向量乘：压缩格式与负载均衡 (对应 Lab 7 / Lab 8: Sparse Matrix Multiply)

// 文件: spmv_dense_baseline.cu
// 用途: 稠密矩阵-向量乘 (GEMV) 基线, 提供"完美规则访存"的性能上限参照
// 编译: nvcc -O3 -arch=sm_80 spmv_dense_baseline.cu -o spmv_dense_baseline
// 运行: ./spmv_dense_baseline 4096 200      (矩阵维度 N, 重复次数 reps)
#include <cstdio>
#include <cstdlib>
#include <cmath>
#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)

#define TPB 256          // 每 block 线程数
#define RPT 4            // 每线程处理的行数 (register coarsening)

// ------------------------------------------------------------------
// (1) 朴素版: A 行主序, 每线程一行, x 直接从全局内存读
//     问题: warp 内相邻线程读的行相隔 N 个 float -> 访存不合并
// ------------------------------------------------------------------
__global__ void gemv_rowmajor(const float * __restrict__ A,
                              const float * __restrict__ x,
                              float * __restrict__ y, int N)
{
    int row = blockIdx.x * blockDim.x + threadIdx.x;
    if (row < N) {
        float dot = 0.0f;
        const float *rp = A + (size_t)row * N;
        for (int j = 0; j < N; ++j) {
            dot = fmaf(rp[j], x[j], dot);
        }
        y[row] = dot;
    }
}

// ------------------------------------------------------------------
// (2) 优化版: A 列主序 (Acol[j*N + i]), 每线程 RPT 行, 可选共享内存缓存 x
//     关键: 固定 j 时, warp 内 32 个 lane 读 Acol[j*N + base + tid] -> 连续 128 B
// ------------------------------------------------------------------
__global__ void gemv_colmajor(const float * __restrict__ Acol,
                              const float * __restrict__ xg,
                              float * __restrict__ y, int N, int cacheX)
{
    extern __shared__ float xs[];
    if (cacheX) {
        for (int j = threadIdx.x; j < N; j += blockDim.x) {
            xs[j] = xg[j];
        }
        __syncthreads();
    }
    const float * __restrict__ x = cacheX ? xs : xg;

    const int bd   = blockDim.x;
    const int base = blockIdx.x * (bd * RPT) + threadIdx.x;

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

    for (int j = 0; j < N; ++j) {
        const float xj = x[j];
        const float *col = Acol + (size_t)j * N + base;
#pragma unroll
        for (int r = 0; r < RPT; ++r) {
            acc[r] = fmaf(col[r * bd], xj, acc[r]);
        }
    }
#pragma unroll
    for (int r = 0; r < RPT; ++r) {
        y[base + r * bd] = acc[r];
    }
}

// ------------------------------------------------------------------
// 计时辅助: 重复 reps 次 kernel, 返回平均每次耗时 (ms)
// ------------------------------------------------------------------
static float time_gemv_rowmajor(const float *d_A, const float *d_x, float *d_y,
                                int N, int reps, float *gflops, float *gbs)
{
    cudaEvent_t s, e;
    CUDA_CHECK(cudaEventCreate(&s));
    CUDA_CHECK(cudaEventCreate(&e));
    gemv_rowmajor<<<(N + TPB - 1) / TPB, TPB>>>(d_A, d_x, d_y, N);   // 预热
    CUDA_CHECK(cudaDeviceSynchronize());
    CUDA_CHECK(cudaEventRecord(s));
    for (int i = 0; i < reps; ++i) {
        gemv_rowmajor<<<(N + TPB - 1) / TPB, TPB>>>(d_A, d_x, d_y, N);
    }
    CUDA_CHECK(cudaEventRecord(e));
    CUDA_CHECK(cudaEventSynchronize(e));
    float ms = 0.0f;
    CUDA_CHECK(cudaEventElapsedTime(&ms, s, e));
    ms /= (float)reps;
    double bytes = 4.0 * (double)N * (double)N + 8.0 * (double)N;
    *gflops = (2.0 * (double)N * (double)N) / (ms * 1.0e6);
    *gbs    = bytes / (ms * 1.0e6);
    CUDA_CHECK(cudaEventDestroy(s));
    CUDA_CHECK(cudaEventDestroy(e));
    return ms;
}

static float time_gemv_colmajor(const float *d_At, const float *d_x, float *d_y,
                                int N, int reps, int cacheX, float *gflops,
                                float *gbs)
{
    const int grid  = N / (TPB * RPT);
    const int smem  = cacheX ? (int)(N * sizeof(float)) : 0;
    cudaEvent_t s, e;
    CUDA_CHECK(cudaEventCreate(&s));
    CUDA_CHECK(cudaEventCreate(&e));
    gemv_colmajor<<<grid, TPB, smem>>>(d_At, d_x, d_y, N, cacheX);   // 预热
    CUDA_CHECK(cudaDeviceSynchronize());
    CUDA_CHECK(cudaEventRecord(s));
    for (int i = 0; i < reps; ++i) {
        gemv_colmajor<<<grid, TPB, smem>>>(d_At, d_x, d_y, N, cacheX);
    }
    CUDA_CHECK(cudaEventRecord(e));
    CUDA_CHECK(cudaEventSynchronize(e));
    float ms = 0.0f;
    CUDA_CHECK(cudaEventElapsedTime(&ms, s, e));
    ms /= (float)reps;
    double bytes = 4.0 * (double)N * (double)N + 8.0 * (double)N;
    *gflops = (2.0 * (double)N * (double)N) / (ms * 1.0e6);
    *gbs    = bytes / (ms * 1.0e6);
    CUDA_CHECK(cudaEventDestroy(s));
    CUDA_CHECK(cudaEventDestroy(e));
    return ms;
}

int main(int argc, char **argv)
{
    const int Nreq = (argc > 1) ? atoi(argv[1]) : 4096;
    const int reps = (argc > 2) ? atoi(argv[2]) : 200;
    // 补齐到 TPB*RPT 的整数倍, 使 kernel 内无需边界判断 (数据用 0 填充)
    const int N = ((Nreq + TPB * RPT - 1) / (TPB * RPT)) * (TPB * RPT);

    const size_t aElems = (size_t)N * (size_t)N;
    printf("N = %d (请求 %d), A = %.1f MB, reps = %d\n",
           N, Nreq, aElems * sizeof(float) / 1.0e6, reps);

    float *h_A   = (float *)malloc(aElems * sizeof(float));
    float *h_At  = (float *)malloc(aElems * sizeof(float));
    float *h_x   = (float *)malloc((size_t)N * sizeof(float));
    float *h_y   = (float *)malloc((size_t)N * sizeof(float));
    float *h_ref = (float *)malloc((size_t)N * sizeof(float));
    if (!h_A || !h_At || !h_x || !h_y || !h_ref) {
        fprintf(stderr, "host malloc failed\n");
        return EXIT_FAILURE;
    }

    srand(12345);
    for (int i = 0; i < N; ++i) h_x[i] = (float)rand() / (float)RAND_MAX - 0.5f;
    for (int i = 0; i < N; ++i) {
        for (int j = 0; j < N; ++j) {
            float v = (float)rand() / (float)RAND_MAX - 0.5f;
            h_A[(size_t)i * N + j]  = v;
            h_At[(size_t)j * N + i] = v;      // 同一份数据的列主序拷贝
        }
    }

    // CPU 参考实现 (double 累加, 减少浮点误差)
    for (int i = 0; i < N; ++i) {
        double dot = 0.0;
        for (int j = 0; j < N; ++j) dot += (double)h_A[(size_t)i * N + j] * h_x[j];
        h_ref[i] = (float)dot;
    }

    float *d_A, *d_At, *d_x, *d_y;
    CUDA_CHECK(cudaMalloc((void **)&d_A,  aElems * sizeof(float)));
    CUDA_CHECK(cudaMalloc((void **)&d_At, aElems * sizeof(float)));
    CUDA_CHECK(cudaMalloc((void **)&d_x,  (size_t)N * sizeof(float)));
    CUDA_CHECK(cudaMalloc((void **)&d_y,  (size_t)N * sizeof(float)));
    CUDA_CHECK(cudaMemcpy(d_A,  h_A,  aElems * sizeof(float), cudaMemcpyHostToDevice));
    CUDA_CHECK(cudaMemcpy(d_At, h_At, aElems * sizeof(float), cudaMemcpyHostToDevice));
    CUDA_CHECK(cudaMemcpy(d_x,  h_x,  (size_t)N * sizeof(float), cudaMemcpyHostToDevice));

    printf("\n%-34s %10s %10s %10s %8s\n",
           "kernel", "time(us)", "GB/s", "GFLOP/s", "verify");
    printf("---------------------------------------------------------------"
           "-------------\n");

    // ---- 版本 1: 行主序朴素版 ----
    float gf = 0.0f, gb = 0.0f;
    float ms = time_gemv_rowmajor(d_A, d_x, d_y, N, reps, &gf, &gb);
    CUDA_CHECK(cudaMemcpy(h_y, d_y, (size_t)N * sizeof(float), cudaMemcpyDeviceToHost));
    int ok = 1;
    for (int i = 0; i < Nreq; ++i) {
        if (fabsf(h_y[i] - h_ref[i]) > 1e-3f * (1.0f + fabsf(h_ref[i]))) { ok = 0; break; }
    }
    printf("%-34s %10.1f %10.1f %10.1f %8s\n",
           "gemv_rowmajor (每线程一行)", ms * 1000.0f, gb, gf, ok ? "PASS" : "FAIL");

    // ---- 版本 2: 列主序 + 寄存器粗化, x 从全局内存读 ----
    ms = time_gemv_colmajor(d_At, d_x, d_y, N, reps, 0, &gf, &gb);
    CUDA_CHECK(cudaMemcpy(h_y, d_y, (size_t)N * sizeof(float), cudaMemcpyDeviceToHost));
    ok = 1;
    for (int i = 0; i < Nreq; ++i) {
        if (fabsf(h_y[i] - h_ref[i]) > 1e-3f * (1.0f + fabsf(h_ref[i]))) { ok = 0; break; }
    }
    printf("%-34s %10.1f %10.1f %10.1f %8s\n",
           "gemv_colmajor (x 读全局)", ms * 1000.0f, gb, gf, ok ? "PASS" : "FAIL");

    // ---- 版本 3: 列主序 + 寄存器粗化 + 共享内存缓存 x ----
    const int smemOK = (N * (int)sizeof(float) <= 48 * 1024) ? 1 : 0;
    ms = time_gemv_colmajor(d_At, d_x, d_y, N, reps, smemOK, &gf, &gb);
    CUDA_CHECK(cudaMemcpy(h_y, d_y, (size_t)N * sizeof(float), cudaMemcpyDeviceToHost));
    ok = 1;
    for (int i = 0; i < Nreq; ++i) {
        if (fabsf(h_y[i] - h_ref[i]) > 1e-3f * (1.0f + fabsf(h_ref[i]))) { ok = 0; break; }
    }
    printf("%-34s %10.1f %10.1f %10.1f %8s\n",
           "gemv_colmajor (x 放共享内存)", ms * 1000.0f, gb, gf, ok ? "PASS" : "FAIL");

    printf("\n参考: A100 峰值 1555 GB/s, 19.5 TFLOPS FP32, 机器平衡点 12.5 FLOP/B\n");
    printf("GEMV 算术强度 = 2*N^2 / (4*N^2 + 8*N) ≈ 0.50 FLOP/B -> Roofline 上限 %.1f GFLOP/s\n",
           1555.0 * 0.5);

    CUDA_CHECK(cudaFree(d_A));
    CUDA_CHECK(cudaFree(d_At));
    CUDA_CHECK(cudaFree(d_x));
    CUDA_CHECK(cudaFree(d_y));
    free(h_A); free(h_At); free(h_x); free(h_y); free(h_ref);
    return EXIT_SUCCESS;
}

// 文件: spmv_csr.cu
// 用途: CSR 格式 SpMV (每线程一行), 含稠密矩阵 -> CSR 的转换函数、CPU 参考实现、
//       计时、有效带宽统计与 warp 负载不均衡分析
// 编译: nvcc -O3 -arch=sm_80 spmv_csr.cu -o spmv_csr
// 运行: ./spmv_csr 1048576 100          (大矩阵行数 N, 重复次数 reps)
#include <cstdio>
#include <cstdlib>
#include <cmath>
#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)

#define TPB 256

// ==================================================================
// 主机端工具 1: 稠密矩阵 -> CSR            (CSR 构建的标准三步法)
//   A: N x N 行主序稠密矩阵 (0 表示稀疏结构的"空位")
//   输出 rowPtr[N+1], colIdx[nnz], vals[nnz]
// ==================================================================
void dense_to_csr(const float *A, int N,
                  int **rowPtr_out, int **colIdx_out, float **vals_out,
                  int *nnz_out)
{
    int nnz = 0;
    for (size_t i = 0; i < (size_t)N * (size_t)N; ++i) {
        if (A[i] != 0.0f) ++nnz;
    }
    const size_t cap = (nnz > 0) ? (size_t)nnz : 1u;

    int   *rowPtr = (int *)malloc(sizeof(int) * (size_t)(N + 1));
    int   *colIdx = (int *)malloc(sizeof(int) * cap);
    float *vals   = (float *)malloc(sizeof(float) * cap);
    if (!rowPtr || !colIdx || !vals) {
        fprintf(stderr, "host malloc failed\n");
        exit(EXIT_FAILURE);
    }

    int k = 0;
    rowPtr[0] = 0;
    for (int i = 0; i < N; ++i) {
        for (int j = 0; j < N; ++j) {
            const float v = A[(size_t)i * (size_t)N + j];
            if (v != 0.0f) {
                colIdx[k] = j;
                vals[k]   = v;
                ++k;
            }
        }
        rowPtr[i + 1] = k;          // 前缀和: 第 i 行的元素区间为 rowPtr[i] 到 rowPtr[i+1] 的左闭右开区间
    }
    *rowPtr_out = rowPtr;
    *colIdx_out = colIdx;
    *vals_out   = vals;
    *nnz_out    = nnz;
}

// ==================================================================
// 主机端工具 2: 生成 N x N 稠密 9 点模板矩阵 (仅用于小规模正确性验证)
//   行 r 对应 2D 网格 (ny x nx) 上的点, 非零列是它自己与 8 个邻居
// ==================================================================
void gen_dense_stencil9(float *A, int N, int nx)
{
    for (size_t i = 0; i < (size_t)N * (size_t)N; ++i) A[i] = 0.0f;
    const int ny = N / nx;
    for (int r = 0; r < N; ++r) {
        const int y = r / nx, xc = r % nx;
        for (int dy = -1; dy <= 1; ++dy) {
            for (int dx = -1; dx <= 1; ++dx) {
                const int yy = y + dy, xx = xc + dx;
                if (yy < 0 || yy >= ny || xx < 0 || xx >= nx) continue;
                const int c = yy * nx + xx;
                A[(size_t)r * (size_t)N + c] =
                    ((c == r) ? 4.0f : 1.0f) + 0.001f * (float)((r * 31 + c * 17) % 13);
            }
        }
    }
}

// ==================================================================
// 主机端工具 3: 直接生成大矩阵的 CSR (不物化 N x N 稠密矩阵)
//   N = 10^6 时稠密矩阵需要 4 TB, 所以大矩阵只能直接构造 CSR
// ==================================================================
void build_stencil9_csr(int numRows, int nx,
                        int **rowPtr_out, int **colIdx_out, float **vals_out,
                        int *nnz_out)
{
    const int ny = numRows / nx;
    int *rowPtr = (int *)malloc(sizeof(int) * (size_t)(numRows + 1));
    if (!rowPtr) { fprintf(stderr, "host malloc failed\n"); exit(EXIT_FAILURE); }

    rowPtr[0] = 0;
    for (int r = 0; r < numRows; ++r) {
        const int y = r / nx, xc = r % nx;
        int cnt = 0;
        for (int dy = -1; dy <= 1; ++dy) {
            for (int dx = -1; dx <= 1; ++dx) {
                const int yy = y + dy, xx = xc + dx;
                if (yy >= 0 && yy < ny && xx >= 0 && xx < nx) ++cnt;
            }
        }
        rowPtr[r + 1] = rowPtr[r] + cnt;
    }
    const int nnz = rowPtr[numRows];
    int   *colIdx = (int *)malloc(sizeof(int) * (size_t)nnz);
    float *vals   = (float *)malloc(sizeof(float) * (size_t)nnz);
    if (!colIdx || !vals) { fprintf(stderr, "host malloc failed\n"); exit(EXIT_FAILURE); }

    int k = 0;
    for (int r = 0; r < numRows; ++r) {
        const int y = r / nx, xc = r % nx;
        for (int dy = -1; dy <= 1; ++dy) {
            for (int dx = -1; dx <= 1; ++dx) {
                const int yy = y + dy, xx = xc + dx;
                if (yy < 0 || yy >= ny || xx < 0 || xx >= nx) continue;
                const int c = yy * nx + xx;
                colIdx[k] = c;
                vals[k]   = ((c == r) ? 4.0f : 1.0f) + 0.001f * (float)((r * 31 + c * 17) % 13);
                ++k;
            }
        }
    }
    *rowPtr_out = rowPtr;
    *colIdx_out = colIdx;
    *vals_out   = vals;
    *nnz_out    = nnz;
}

// ==================================================================
// 主机端 CPU 参考实现 (直接按 CSR 语义计算)
// ==================================================================
void spmv_csr_cpu(int numRows, const int *rowPtr, const int *colIdx,
                  const float *vals, const float *x, float *y)
{
    for (int r = 0; r < numRows; ++r) {
        double dot = 0.0;
        for (int e = rowPtr[r]; e < rowPtr[r + 1]; ++e) {
            dot += (double)vals[e] * (double)x[colIdx[e]];
        }
        y[r] = (float)dot;
    }
}

// ==================================================================
// GPU kernel: 讲义 page 15 的 SpMV_CSR, 补上 __ldg 与 fmaf
// ==================================================================
__global__ void spmv_csr(int num_rows,
                         const float * __restrict__ data,
                         const int   * __restrict__ col_index,
                         const int   * __restrict__ row_ptr,
                         const float * __restrict__ x,
                         float       * __restrict__ y)
{
    int row = blockIdx.x * blockDim.x + threadIdx.x;
    if (row < num_rows) {
        float dot = 0.0f;
        int row_start = __ldg(&row_ptr[row]);
        int row_end   = __ldg(&row_ptr[row + 1]);
        for (int elem = row_start; elem < row_end; elem++) {
            dot = fmaf(__ldg(&data[elem]), __ldg(&x[__ldg(&col_index[elem])]), dot);
        }
        y[row] = dot;
    }
}

int main(int argc, char **argv)
{
    // ==============================================================
    // 第 1 部分: 稠密 -> CSR 转换与正确性验证 (小矩阵, 稠密可以物化)
    // ==============================================================
    const int Ns = 1024, nxs = 32;
    printf("=== 第 1 部分: dense_to_csr 与 kernel 正确性验证 (N = %d) ===\n", Ns);

    float *h_Ad = (float *)calloc((size_t)Ns * (size_t)Ns, sizeof(float));
    if (!h_Ad) { fprintf(stderr, "host malloc failed\n"); return EXIT_FAILURE; }
    gen_dense_stencil9(h_Ad, Ns, nxs);

    int   *s_rowPtr = NULL, *s_colIdx = NULL, s_nnz = 0;
    float *s_vals   = NULL;
    dense_to_csr(h_Ad, Ns, &s_rowPtr, &s_colIdx, &s_vals, &s_nnz);

    const double denseBytes = (double)Ns * Ns * sizeof(float);
    const double csrBytes   = 8.0 * s_nnz + 4.0 * (Ns + 1);
    printf("稠密存储 = %.2f MB,  CSR 存储 = 8*nnz + 4*(N+1) = %.3f MB,  nnz = %d\n",
           denseBytes / 1.0e6, csrBytes / 1.0e6, s_nnz);
    printf("CSR/稠密 = %.2f%%  (矩阵越稀疏, 压缩收益越大)\n",
           100.0 * csrBytes / denseBytes);

    float *s_x = (float *)malloc(sizeof(float) * (size_t)Ns);
    float *s_y = (float *)malloc(sizeof(float) * (size_t)Ns);
    float *s_ref = (float *)malloc(sizeof(float) * (size_t)Ns);
    if (!s_x || !s_y || !s_ref) { fprintf(stderr, "host malloc failed\n"); return EXIT_FAILURE; }
    for (int i = 0; i < Ns; ++i) s_x[i] = (float)((i % 17) - 8) * 0.25f;

    spmv_csr_cpu(Ns, s_rowPtr, s_colIdx, s_vals, s_x, s_ref);

    int   *d_srow = NULL, *d_scol = NULL;
    float *d_sval = NULL, *d_sx = NULL, *d_sy = NULL;
    CUDA_CHECK(cudaMalloc((void **)&d_srow, sizeof(int) * (size_t)(Ns + 1)));
    CUDA_CHECK(cudaMalloc((void **)&d_scol, sizeof(int) * (size_t)s_nnz));
    CUDA_CHECK(cudaMalloc((void **)&d_sval, sizeof(float) * (size_t)s_nnz));
    CUDA_CHECK(cudaMalloc((void **)&d_sx, sizeof(float) * (size_t)Ns));
    CUDA_CHECK(cudaMalloc((void **)&d_sy, sizeof(float) * (size_t)Ns));
    CUDA_CHECK(cudaMemcpy(d_srow, s_rowPtr, sizeof(int) * (size_t)(Ns + 1), cudaMemcpyHostToDevice));
    CUDA_CHECK(cudaMemcpy(d_scol, s_colIdx, sizeof(int) * (size_t)s_nnz, cudaMemcpyHostToDevice));
    CUDA_CHECK(cudaMemcpy(d_sval, s_vals, sizeof(float) * (size_t)s_nnz, cudaMemcpyHostToDevice));
    CUDA_CHECK(cudaMemcpy(d_sx, s_x, sizeof(float) * (size_t)Ns, cudaMemcpyHostToDevice));

    spmv_csr<<<(Ns + TPB - 1) / TPB, TPB>>>(Ns, d_sval, d_scol, d_srow, d_sx, d_sy);
    CUDA_CHECK(cudaGetLastError());
    CUDA_CHECK(cudaDeviceSynchronize());
    CUDA_CHECK(cudaMemcpy(s_y, d_sy, sizeof(float) * (size_t)Ns, cudaMemcpyDeviceToHost));

    int ok = 1;
    for (int i = 0; i < Ns; ++i) {
        if (fabsf(s_y[i] - s_ref[i]) > 1e-4f * (1.0f + fabsf(s_ref[i]))) { ok = 0; break; }
    }
    printf("CSR SpMV 验证: %s\n\n", ok ? "PASS" : "FAIL");

    CUDA_CHECK(cudaFree(d_srow)); CUDA_CHECK(cudaFree(d_scol));
    CUDA_CHECK(cudaFree(d_sval)); CUDA_CHECK(cudaFree(d_sx));
    CUDA_CHECK(cudaFree(d_sy));
    free(h_Ad); free(s_rowPtr); free(s_colIdx); free(s_vals);
    free(s_x); free(s_y); free(s_ref);

    // ==============================================================
    // 第 2 部分: 大矩阵性能测试 (直接构造 CSR, 不物化稠密矩阵)
    // ==============================================================
    const int N    = (argc > 1) ? atoi(argv[1]) : 1048576;
    const int reps = (argc > 2) ? atoi(argv[2]) : 100;

    int   *rowPtr = NULL, *colIdx = NULL, nnz = 0;
    float *vals   = NULL;
    build_stencil9_csr(N, 1024, &rowPtr, &colIdx, &vals, &nnz);

    double maxLen = 0.0, sumLen = 0.0;
    for (int r = 0; r < N; ++r) {
        const double L = rowPtr[r + 1] - rowPtr[r];
        if (L > maxLen) maxLen = L;
        sumLen += L;
    }
    const double avgLen = sumLen / (double)N;

    const double bytes = 4.0 * nnz              // values
                       + 4.0 * nnz              // colIdx
                       + 4.0 * (double)(N + 1)  // rowPtr
                       + 4.0 * (double)N        // x
                       + 4.0 * (double)N;       // y
    printf("=== 第 2 部分: CSR SpMV 性能 (N = %d, nnz = %d) ===\n", N, nnz);
    printf("行长: 平均 %.2f, 最长 %.0f;  必须搬运的字节 = %.1f MB\n",
           avgLen, maxLen, bytes / 1.0e6);
    printf("  其中 values %.1f MB + colIdx %.1f MB + rowPtr %.1f MB + x %.1f MB + y %.1f MB\n",
           4.0 * nnz / 1.0e6, 4.0 * nnz / 1.0e6, 4.0 * (N + 1) / 1.0e6,
           4.0 * N / 1.0e6, 4.0 * N / 1.0e6);

    float *x = (float *)malloc(sizeof(float) * (size_t)N);
    float *ref = (float *)malloc(sizeof(float) * (size_t)N);
    float *y = (float *)malloc(sizeof(float) * (size_t)N);
    if (!x || !ref || !y) { fprintf(stderr, "host malloc failed\n"); return EXIT_FAILURE; }
    for (int i = 0; i < N; ++i) x[i] = (float)((i % 17) - 8) * 0.25f;

    spmv_csr_cpu(N, rowPtr, colIdx, vals, x, ref);

    int   *d_row = NULL, *d_col = NULL;
    float *d_val = NULL, *d_x = NULL, *d_y = NULL;
    CUDA_CHECK(cudaMalloc((void **)&d_row, sizeof(int) * (size_t)(N + 1)));
    CUDA_CHECK(cudaMalloc((void **)&d_col, sizeof(int) * (size_t)nnz));
    CUDA_CHECK(cudaMalloc((void **)&d_val, sizeof(float) * (size_t)nnz));
    CUDA_CHECK(cudaMalloc((void **)&d_x, sizeof(float) * (size_t)N));
    CUDA_CHECK(cudaMalloc((void **)&d_y, sizeof(float) * (size_t)N));
    CUDA_CHECK(cudaMemcpy(d_row, rowPtr, sizeof(int) * (size_t)(N + 1), cudaMemcpyHostToDevice));
    CUDA_CHECK(cudaMemcpy(d_col, colIdx, sizeof(int) * (size_t)nnz, cudaMemcpyHostToDevice));
    CUDA_CHECK(cudaMemcpy(d_val, vals, sizeof(float) * (size_t)nnz, cudaMemcpyHostToDevice));
    CUDA_CHECK(cudaMemcpy(d_x, x, sizeof(float) * (size_t)N, cudaMemcpyHostToDevice));

    const int grid = (N + TPB - 1) / TPB;

    spmv_csr<<<grid, TPB>>>(N, d_val, d_col, d_row, d_x, d_y);
    CUDA_CHECK(cudaGetLastError());
    CUDA_CHECK(cudaDeviceSynchronize());
    CUDA_CHECK(cudaMemcpy(y, d_y, sizeof(float) * (size_t)N, cudaMemcpyDeviceToHost));
    ok = 1;
    for (int i = 0; i < N; ++i) {
        if (fabsf(y[i] - ref[i]) > 1e-4f * (1.0f + fabsf(ref[i]))) { ok = 0; break; }
    }
    printf("大矩阵验证: %s\n", ok ? "PASS" : "FAIL");

    cudaEvent_t s, e;
    CUDA_CHECK(cudaEventCreate(&s));
    CUDA_CHECK(cudaEventCreate(&e));
    spmv_csr<<<grid, TPB>>>(N, d_val, d_col, d_row, d_x, d_y);
    CUDA_CHECK(cudaDeviceSynchronize());
    CUDA_CHECK(cudaEventRecord(s));
    for (int i = 0; i < reps; ++i) {
        spmv_csr<<<grid, TPB>>>(N, d_val, d_col, d_row, d_x, d_y);
    }
    CUDA_CHECK(cudaEventRecord(e));
    CUDA_CHECK(cudaEventSynchronize(e));
    float ms = 0.0f;
    CUDA_CHECK(cudaEventElapsedTime(&ms, s, e));
    ms /= (float)reps;

    const double gflops = (2.0 * nnz) / (ms * 1.0e6);
    const double gbs    = bytes / (ms * 1.0e6);
    printf("\nCSR SpMV (每线程一行): %.1f us, 有效带宽 %.1f GB/s (峰值 %.1f GB/s 的 %.1f%%), "
           "%.1f GFLOP/s (峰值的 %.2f%%)\n",
           ms * 1000.0f, gbs, 1555.0, 100.0 * gbs / 1555.0, gflops,
           100.0 * gflops / 19500.0);
    printf("说明: 有效带宽 = 必须搬运的字节 / 时间 (含 colIdx 与 rowPtr 的开销);\n");
    printf("      gather 造成的额外 sector 放大不计入分子, 所以它是\"有用字节口径\"的带宽。\n");

    CUDA_CHECK(cudaFree(d_row)); CUDA_CHECK(cudaFree(d_col)); CUDA_CHECK(cudaFree(d_val));
    CUDA_CHECK(cudaFree(d_x));   CUDA_CHECK(cudaFree(d_y));
    CUDA_CHECK(cudaEventDestroy(s)); CUDA_CHECK(cudaEventDestroy(e));
    free(rowPtr); free(colIdx); free(vals); free(x); free(ref); free(y);
    return EXIT_SUCCESS;
}

// 文件: spmv_ell.cu
// 用途: ELL(PACK) 格式 SpMV: CSR -> ELL 转换、基础 ELL kernel、float4 向量化 ELL kernel,
//       并与 CSR 版本在同一矩阵上实测对比; 含"离群长行导致 padding 爆炸"的实验
// 编译: nvcc -O3 -arch=sm_80 spmv_ell.cu -o spmv_ell
// 运行: ./spmv_ell 1048576 100       (大矩阵行数 N, 重复次数 reps)
#include <cstdio>
#include <cstdlib>
#include <cmath>
#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)

#define TPB 256

// ==================================================================
// 主机端: 构造 CSR。hubStride > 0 时, 每隔 hubStride 行额外加入 hubExtra 个列,
//          用来模拟"枢纽节点 / 离群长行"(图矩阵中常见的度分布长尾)
// ==================================================================
void build_stencil9_csr(int numRows, int nx, int hubStride, int hubExtra,
                        int **rowPtr_out, int **colIdx_out, float **vals_out,
                        int *nnz_out)
{
    const int ny = numRows / nx;
    int *rowPtr = (int *)malloc(sizeof(int) * (size_t)(numRows + 1));
    if (!rowPtr) { fprintf(stderr, "host malloc failed\n"); exit(EXIT_FAILURE); }

    rowPtr[0] = 0;
    for (int r = 0; r < numRows; ++r) {
        const int y = r / nx, xc = r % nx;
        int cnt = 0;
        for (int dy = -1; dy <= 1; ++dy) {
            for (int dx = -1; dx <= 1; ++dx) {
                const int yy = y + dy, xx = xc + dx;
                if (yy >= 0 && yy < ny && xx >= 0 && xx < nx) ++cnt;
            }
        }
        if (hubStride > 0 && (r % hubStride) == 0) cnt += hubExtra;
        rowPtr[r + 1] = rowPtr[r] + cnt;
    }
    const int nnz = rowPtr[numRows];
    int   *colIdx = (int *)malloc(sizeof(int) * (size_t)nnz);
    float *vals   = (float *)malloc(sizeof(float) * (size_t)nnz);
    if (!colIdx || !vals) { fprintf(stderr, "host malloc failed\n"); exit(EXIT_FAILURE); }

    int k = 0;
    for (int r = 0; r < numRows; ++r) {
        const int y = r / nx, xc = r % nx;
        for (int dy = -1; dy <= 1; ++dy) {
            for (int dx = -1; dx <= 1; ++dx) {
                const int yy = y + dy, xx = xc + dx;
                if (yy < 0 || yy >= ny || xx < 0 || xx >= nx) continue;
                const int c = yy * nx + xx;
                colIdx[k] = c;
                vals[k]   = ((c == r) ? 4.0f : 1.0f) + 0.001f * (float)((r * 31 + c * 17) % 13);
                ++k;
            }
        }
        if (hubStride > 0 && (r % hubStride) == 0) {
            for (int e = 0; e < hubExtra; ++e) {
                const int c = (int)(((long long)r * 7919LL + (long long)e * 104729LL) % numRows);
                colIdx[k] = c;
                vals[k]   = 0.5f + 0.001f * (float)((r + e * 7) % 11);
                ++k;
            }
        }
    }
    *rowPtr_out = rowPtr;
    *colIdx_out = colIdx;
    *vals_out   = vals;
    *nnz_out    = nnz;
}

// ==================================================================
// 主机端: CSR -> ELL (列主序布局: data[i*numRows + row], 不足处补 0)
//   E = 最长行的长度; 返回 E 与实际分配的字节数统计
// ==================================================================
void csr_to_ell(int numRows, const int *rowPtr, const int *colIdx, const float *vals,
                int *numElem_out, int **ellCol_out, float **ellVal_out,
                double *padRatio_out)
{
    int E = 0;
    for (int r = 0; r < numRows; ++r) {
        const int L = rowPtr[r + 1] - rowPtr[r];
        if (L > E) E = L;
    }
    const size_t total = (size_t)numRows * (size_t)E;
    int   *ellCol = (int *)malloc(sizeof(int) * total);
    float *ellVal = (float *)malloc(sizeof(float) * total);
    if (!ellCol || !ellVal) { fprintf(stderr, "host malloc failed\n"); exit(EXIT_FAILURE); }

    for (size_t i = 0; i < total; ++i) {       // 先整体填充 padding: 值 0, 列号 0
        ellCol[i] = 0;
        ellVal[i] = 0.0f;
    }
    for (int r = 0; r < numRows; ++r) {
        const int start = rowPtr[r];
        const int len   = rowPtr[r + 1] - rowPtr[r];
        for (int i = 0; i < len; ++i) {
            ellVal[(size_t)i * (size_t)numRows + (size_t)r] = vals[start + i];
            ellCol[(size_t)i * (size_t)numRows + (size_t)r] = colIdx[start + i];
        }
    }
    *numElem_out  = E;
    *ellCol_out   = ellCol;
    *ellVal_out   = ellVal;
    *padRatio_out = (double)total / (double)rowPtr[numRows];
}

// ==================================================================
// 主机端 CPU 参考实现 (按 CSR 语义计算)
// ==================================================================
void spmv_cpu(int numRows, const int *rowPtr, const int *colIdx,
              const float *vals, const float *x, float *y)
{
    for (int r = 0; r < numRows; ++r) {
        double dot = 0.0;
        for (int e = rowPtr[r]; e < rowPtr[r + 1]; ++e) {
            dot += (double)vals[e] * (double)x[colIdx[e]];
        }
        y[r] = (float)dot;
    }
}

// ==================================================================
// kernel 1: 讲义 page 20 的 SpMV_ELL (基础版, 标量 4 B 访问)
// ==================================================================
__global__ void spmv_ell(int num_rows,
                         const float * __restrict__ data,
                         const int   * __restrict__ col_index,
                         int num_elem,
                         const float * __restrict__ x,
                         float       * __restrict__ y)
{
    int row = blockIdx.x * blockDim.x + threadIdx.x;
    if (row < num_rows) {
        float dot = 0.0f;
        for (int i = 0; i < num_elem; i++) {
            const int k = i * num_rows + row;
            dot = fmaf(data[k], __ldg(&x[__ldg(&col_index[k])]), dot);
        }
        y[row] = dot;
    }
}

// ==================================================================
// kernel 2: float4 向量化 ELL, 每线程处理 4 个连续行
//   data 被当作 float4*, 第 i 次迭代的 float4 覆盖 row0..row0+3
//   要求: num_rows % 4 == 0
// ==================================================================
__global__ void spmv_ell_vec4(int num_rows,
                              const float4 * __restrict__ data4,
                              const int4   * __restrict__ col4,
                              int num_elem,
                              const float * __restrict__ x,
                              float       * __restrict__ y)
{
    const int half = num_rows >> 2;
    const int g = blockIdx.x * blockDim.x + threadIdx.x;
    if (g < half) {
        float4 acc = make_float4(0.0f, 0.0f, 0.0f, 0.0f);
        for (int i = 0; i < num_elem; ++i) {
            const float4 d = data4[(size_t)i * (size_t)half + (size_t)g];
            const int4   c = col4[(size_t)i * (size_t)half + (size_t)g];
            acc.x = fmaf(d.x, __ldg(&x[c.x]), acc.x);
            acc.y = fmaf(d.y, __ldg(&x[c.y]), acc.y);
            acc.z = fmaf(d.z, __ldg(&x[c.z]), acc.z);
            acc.w = fmaf(d.w, __ldg(&x[c.w]), acc.w);
        }
        const int row0 = g * 4;
        y[row0 + 0] = acc.x;
        y[row0 + 1] = acc.y;
        y[row0 + 2] = acc.z;
        y[row0 + 3] = acc.w;
    }
}

// ==================================================================
// kernel 0: CSR 版本 (作为对比基准)
// ==================================================================
__global__ void spmv_csr(int num_rows,
                         const float * __restrict__ data,
                         const int   * __restrict__ col_index,
                         const int   * __restrict__ row_ptr,
                         const float * __restrict__ x,
                         float       * __restrict__ y)
{
    int row = blockIdx.x * blockDim.x + threadIdx.x;
    if (row < num_rows) {
        float dot = 0.0f;
        int row_start = __ldg(&row_ptr[row]);
        int row_end   = __ldg(&row_ptr[row + 1]);
        for (int elem = row_start; elem < row_end; elem++) {
            dot = fmaf(__ldg(&data[elem]), __ldg(&x[__ldg(&col_index[elem])]), dot);
        }
        y[row] = dot;
    }
}

// ==================================================================
// 单个测试用例: 构造 CSR -> 转换 ELL -> 三个 kernel 各自验证与计时
// ==================================================================
static void run_case(const char *name, int N, int nx, int hubStride, int hubExtra,
                     int reps)
{
    int   *rowPtr = NULL, *colIdx = NULL, nnz = 0;
    float *vals   = NULL;
    build_stencil9_csr(N, nx, hubStride, hubExtra, &rowPtr, &colIdx, &vals, &nnz);

    double sumLen = 0.0;
    int    maxLen = 0;
    for (int r = 0; r < N; ++r) {
        const int L = rowPtr[r + 1] - rowPtr[r];
        sumLen += L;
        if (L > maxLen) maxLen = L;
    }
    const double avgLen = sumLen / (double)N;

    int   E = 0, *ellCol = NULL;
    float *ellVal = NULL;
    double padRatio = 0.0;
    csr_to_ell(N, rowPtr, colIdx, vals, &E, &ellCol, &ellVal, &padRatio);

    const double csrBytes = 8.0 * nnz + 4.0 * (double)(N + 1) + 8.0 * (double)N;
    const double ellBytes = 8.0 * (double)N * (double)E + 8.0 * (double)N;

    printf("\n=== %s ===\n", name);
    printf("N = %d, nnz = %d, 平均行长 = %.2f, 最长行 E = %d\n", N, nnz, avgLen, E);
    printf("CSR 必须搬运 %.1f MB; ELL 必须搬运 %.1f MB (padding 放大 %.2fx)\n",
           csrBytes / 1.0e6, ellBytes / 1.0e6, padRatio);
    if (E > 4 * avgLen) {
        printf("警告: E / 平均行长 = %.1f > 4, padding 爆炸, 应该改用 HYB(ELL + COO)\n",
               (double)E / avgLen);
    }

    float *x   = (float *)malloc(sizeof(float) * (size_t)N);
    float *ref = (float *)malloc(sizeof(float) * (size_t)N);
    float *y   = (float *)malloc(sizeof(float) * (size_t)N);
    if (!x || !ref || !y) { fprintf(stderr, "host malloc failed\n"); exit(EXIT_FAILURE); }
    for (int i = 0; i < N; ++i) x[i] = (float)((i % 17) - 8) * 0.25f;
    spmv_cpu(N, rowPtr, colIdx, vals, x, ref);

    int   *d_row = NULL, *d_col = NULL, *d_ellCol = NULL;
    float *d_val = NULL, *d_x = NULL, *d_y = NULL, *d_ellVal = NULL;
    CUDA_CHECK(cudaMalloc((void **)&d_row, sizeof(int) * (size_t)(N + 1)));
    CUDA_CHECK(cudaMalloc((void **)&d_col, sizeof(int) * (size_t)nnz));
    CUDA_CHECK(cudaMalloc((void **)&d_val, sizeof(float) * (size_t)nnz));
    CUDA_CHECK(cudaMalloc((void **)&d_ellCol, sizeof(int) * (size_t)N * (size_t)E));
    CUDA_CHECK(cudaMalloc((void **)&d_ellVal, sizeof(float) * (size_t)N * (size_t)E));
    CUDA_CHECK(cudaMalloc((void **)&d_x, sizeof(float) * (size_t)N));
    CUDA_CHECK(cudaMalloc((void **)&d_y, sizeof(float) * (size_t)N));
    CUDA_CHECK(cudaMemcpy(d_row, rowPtr, sizeof(int) * (size_t)(N + 1), cudaMemcpyHostToDevice));
    CUDA_CHECK(cudaMemcpy(d_col, colIdx, sizeof(int) * (size_t)nnz, cudaMemcpyHostToDevice));
    CUDA_CHECK(cudaMemcpy(d_val, vals, sizeof(float) * (size_t)nnz, cudaMemcpyHostToDevice));
    CUDA_CHECK(cudaMemcpy(d_ellCol, ellCol, sizeof(int) * (size_t)N * (size_t)E, cudaMemcpyHostToDevice));
    CUDA_CHECK(cudaMemcpy(d_ellVal, ellVal, sizeof(float) * (size_t)N * (size_t)E, cudaMemcpyHostToDevice));
    CUDA_CHECK(cudaMemcpy(d_x, x, sizeof(float) * (size_t)N, cudaMemcpyHostToDevice));

    cudaEvent_t ev0, ev1;
    CUDA_CHECK(cudaEventCreate(&ev0));
    CUDA_CHECK(cudaEventCreate(&ev1));

    // ---- 三个 kernel 各跑一次并验证 ----
    spmv_csr<<<(N + TPB - 1) / TPB, TPB>>>(N, d_val, d_col, d_row, d_x, d_y);
    CUDA_CHECK(cudaDeviceSynchronize());
    CUDA_CHECK(cudaMemcpy(y, d_y, sizeof(float) * (size_t)N, cudaMemcpyDeviceToHost));
    int okCsr = 1;
    for (int i = 0; i < N; ++i)
        if (fabsf(y[i] - ref[i]) > 1e-4f * (1.0f + fabsf(ref[i]))) { okCsr = 0; break; }

    spmv_ell<<<(N + TPB - 1) / TPB, TPB>>>(N, d_ellVal, d_ellCol, E, d_x, d_y);
    CUDA_CHECK(cudaDeviceSynchronize());
    CUDA_CHECK(cudaMemcpy(y, d_y, sizeof(float) * (size_t)N, cudaMemcpyDeviceToHost));
    int okEll = 1;
    for (int i = 0; i < N; ++i)
        if (fabsf(y[i] - ref[i]) > 1e-4f * (1.0f + fabsf(ref[i]))) { okEll = 0; break; }

    int okEll4 = -1;                 // -1 表示未测试
    const int half = N >> 2;
    if ((N % 4) == 0) {
        spmv_ell_vec4<<<(half + TPB - 1) / TPB, TPB>>>(
            N, (const float4 *)d_ellVal, (const int4 *)d_ellCol, E, d_x, d_y);
        CUDA_CHECK(cudaDeviceSynchronize());
        CUDA_CHECK(cudaMemcpy(y, d_y, sizeof(float) * (size_t)N, cudaMemcpyDeviceToHost));
        okEll4 = 1;
        for (int i = 0; i < N; ++i)
            if (fabsf(y[i] - ref[i]) > 1e-4f * (1.0f + fabsf(ref[i]))) { okEll4 = 0; break; }
    }

    // ---- 计时 ----
    float ms;
    printf("%-22s %10s %10s %11s %9s\n", "kernel", "time(us)", "GB/s", "GFLOP/s", "verify");

    spmv_csr<<<(N + TPB - 1) / TPB, TPB>>>(N, d_val, d_col, d_row, d_x, d_y);
    CUDA_CHECK(cudaDeviceSynchronize());
    CUDA_CHECK(cudaEventRecord(ev0));
    for (int i = 0; i < reps; ++i)
        spmv_csr<<<(N + TPB - 1) / TPB, TPB>>>(N, d_val, d_col, d_row, d_x, d_y);
    CUDA_CHECK(cudaEventRecord(ev1));
    CUDA_CHECK(cudaEventSynchronize(ev1));
    CUDA_CHECK(cudaEventElapsedTime(&ms, ev0, ev1));
    ms /= (float)reps;
    printf("%-22s %10.1f %10.1f %11.1f %9s\n", "CSR (每线程一行)", ms * 1000.0f,
           csrBytes / (ms * 1.0e6), 2.0 * nnz / (ms * 1.0e6), okCsr ? "PASS" : "FAIL");

    spmv_ell<<<(N + TPB - 1) / TPB, TPB>>>(N, d_ellVal, d_ellCol, E, d_x, d_y);
    CUDA_CHECK(cudaDeviceSynchronize());
    CUDA_CHECK(cudaEventRecord(ev0));
    for (int i = 0; i < reps; ++i)
        spmv_ell<<<(N + TPB - 1) / TPB, TPB>>>(N, d_ellVal, d_ellCol, E, d_x, d_y);
    CUDA_CHECK(cudaEventRecord(ev1));
    CUDA_CHECK(cudaEventSynchronize(ev1));
    CUDA_CHECK(cudaEventElapsedTime(&ms, ev0, ev1));
    ms /= (float)reps;
    printf("%-22s %10.1f %10.1f %11.1f %9s\n", "ELL (标量)", ms * 1000.0f,
           ellBytes / (ms * 1.0e6), 2.0 * nnz / (ms * 1.0e6), okEll ? "PASS" : "FAIL");

    if (okEll4 >= 0) {
        spmv_ell_vec4<<<(half + TPB - 1) / TPB, TPB>>>(
            N, (const float4 *)d_ellVal, (const int4 *)d_ellCol, E, d_x, d_y);
        CUDA_CHECK(cudaDeviceSynchronize());
        CUDA_CHECK(cudaEventRecord(ev0));
        for (int i = 0; i < reps; ++i)
            spmv_ell_vec4<<<(half + TPB - 1) / TPB, TPB>>>(
                N, (const float4 *)d_ellVal, (const int4 *)d_ellCol, E, d_x, d_y);
        CUDA_CHECK(cudaEventRecord(ev1));
        CUDA_CHECK(cudaEventSynchronize(ev1));
        CUDA_CHECK(cudaEventElapsedTime(&ms, ev0, ev1));
        ms /= (float)reps;
        printf("%-22s %10.1f %10.1f %11.1f %9s\n", "ELL + float4", ms * 1000.0f,
               ellBytes / (ms * 1.0e6), 2.0 * nnz / (ms * 1.0e6),
               okEll4 ? "PASS" : "FAIL");
    }

    CUDA_CHECK(cudaFree(d_row));    CUDA_CHECK(cudaFree(d_col));
    CUDA_CHECK(cudaFree(d_val));    CUDA_CHECK(cudaFree(d_ellCol));
    CUDA_CHECK(cudaFree(d_ellVal)); CUDA_CHECK(cudaFree(d_x));
    CUDA_CHECK(cudaFree(d_y));
    CUDA_CHECK(cudaEventDestroy(ev0)); CUDA_CHECK(cudaEventDestroy(ev1));
    free(rowPtr); free(colIdx); free(vals);
    free(ellCol); free(ellVal); free(x); free(ref); free(y);
}

int main(int argc, char **argv)
{
    const int N    = (argc > 1) ? atoi(argv[1]) : 1048576;
    const int reps = (argc > 2) ? atoi(argv[2]) : 100;

    printf("A100 基准: 1555 GB/s, 19.5 TFLOPS FP32, L2 40 MB, 108 SM, 164 KB smem/SM\n");

    // 矩阵 A: 纯 9 点模板 (行长 4/6/9, 方差极小) -> ELL 的最佳场景
    run_case("矩阵 A: 9 点模板 (行长方差极小)", N, 1024, 0, 0, reps);

    // 矩阵 B: 9 点模板 + 每 1024 行注入 40 个额外非零元 (离群长行) -> padding 爆炸
    run_case("矩阵 B: 9 点模板 + 离群长行 (每 1024 行 +40 元)", N, 1024, 1024, 40, reps);

    return EXIT_SUCCESS;
}