Lecture 8: 并行模式六 —— 分块卷积：常量内存与边界处理 (对应 Lab 6 / Lab 4: 3D Convolution)

/* (a) 补零: 越界返回 0 —— 课程 Lab 的默认约定 */
__device__ __forceinline__ float get_zero(const float *N, int idx, int len) {
    return (idx >= 0 && idx < len) ? N[idx] : 0.0f;
}
/* (b) 夹紧: 越界贴到最近的合法下标 */
__device__ __forceinline__ float get_clamp(const float *N, int idx, int len) {
    int k = idx < 0 ? 0 : (idx >= len ? len - 1 : idx);
    return N[k];
}
/* (c) 环绕: 用取模把下标折回数组内 (len 必须是编译期可知或用 % 运算) */
__device__ __forceinline__ float get_wrap(const float *N, int idx, int len) {
    int k = idx % len;
    if (k < 0) k += len;
    return N[k];
}

for (int j = 0; j < MASK_WIDTH; ++j) {
    int idx = N_start_point + j;
    if (idx >= 0 && idx < Width)          /* 这一句就是发散源 */
        Pvalue += N[idx] * Mc[j];
}

float v = 0.0f;
if (row_i >= 0 && row_i < Height && col_i >= 0 && col_i < Width)
    v = N[row_i * Width + col_i];
Ns[ty][tx] = v;                  /* 所有线程都执行这一句, 分支外, 安全 */
__syncthreads();                 /* __syncthreads 必须在分支之外! */

/* 方式一: 一维展平 (示例 3 采用, 简单直观)
   idx = x + y*W + z*W*H;  cudaMalloc(&d_N, W*H*D*sizeof(float));
   一次 cudaMemcpy 传输 W*H*D 个 float。 */

/* 方式二: cudaMalloc3D + cudaPitchedPtr (对行首对齐更友好) */
cudaExtent extent = make_cudaExtent(W * sizeof(float), H, D);
cudaPitchedPtr d_Np;                       /* 内含 ptr / pitch / xsize / ysize */
CUDA_CHECK(cudaMalloc3D(&d_Np, extent));  /* 硬件会把 pitch 对齐到 512 字节 */
cudaMemcpy3DParms p = {0};
p.srcPtr = make_cudaPitchedPtr((void *)h_N, W * sizeof(float), W, H);
p.dstPtr = d_Np;
p.extent = extent;
p.kind   = cudaMemcpyHostToDevice;
CUDA_CHECK(cudaMemcpy3D(&p));              /* 一次调用完成 W*H*D 的 H2D 拷贝 */
/* kernel 中访问 (x, y, z): 先用 int Wp = d_Np.pitch / sizeof(float); 得到行长,
   然后 index = z * Wp * H + y * Wp + x。 pitch >= W*4, 多出的填充字节不可读。 */

// 文件: conv2d_naive_constant.cu
// 编译: nvcc -O3 -arch=sm_80 conv2d_naive_constant.cu -o conv2d_naive
// 运行: ./conv2d_naive             # 默认 2048x2048 输入, 5x5 mask
//       ./conv2d_naive 1024 1024   # 指定输入尺寸
//
// 功能: 朴素(未分块) 2D 卷积, 每个线程计算一个输出元素。
//       - 5x5 卷积核放在 __constant__ 常量内存 Mc 中, 用 cudaMemcpyToSymbol 初始化
//       - 输出尺寸与输入相同(centered 形式), 越界的 ghost 元素按 0 处理 (zero padding)
//       - 打印耗时 / GB/s / GFLOP/s, 并与 CPU 参考实现逐元素对比

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

#define MASK_WIDTH 5
#define MASK_RADIUS (MASK_WIDTH / 2)

/* 常量内存: 文件作用域声明, 所有 kernel 可见, 上限 64 KB */
static __constant__ float Mc[MASK_WIDTH * MASK_WIDTH];

#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)

/* ------------------------- GPU: 朴素 2D 卷积 ------------------------- */
__global__ void conv2d_naive_kernel(const float * __restrict__ N,
                                    float * __restrict__ P,
                                    int Width, int Height)
{
    int col = blockIdx.x * blockDim.x + threadIdx.x;
    int row = blockIdx.y * blockDim.y + threadIdx.y;
    if (row >= Height || col >= Width) return;      /* 输出越界线程直接退出 */

    float Pvalue = 0.0f;
    #pragma unroll
    for (int i = 0; i < MASK_WIDTH; ++i) {
        int nrow = row - MASK_RADIUS + i;
        if (nrow < 0 || nrow >= Height) continue;   /* ghost 行 -> 视作 0 */
        #pragma unroll
        for (int j = 0; j < MASK_WIDTH; ++j) {
            int ncol = col - MASK_RADIUS + j;
            if (ncol < 0 || ncol >= Width) continue;/* ghost 列 -> 视作 0 */
            Pvalue += N[nrow * Width + ncol] * Mc[i * MASK_WIDTH + j];
        }
    }
    P[row * Width + col] = Pvalue;
}

/* 主机端 mask 副本, 供 CPU 参考实现使用 (与 device 端 Mc 内容一致) */
static float h_Mck[MASK_WIDTH * MASK_WIDTH];

/* ------------------------- CPU 参考实现 ------------------------- */
static void conv2d_cpu_reference(const float *N, float *P, int Width, int Height)
{
    for (int row = 0; row < Height; ++row) {
        for (int col = 0; col < Width; ++col) {
            float sum = 0.0f;
            for (int i = 0; i < MASK_WIDTH; ++i) {
                int nrow = row - MASK_RADIUS + i;
                if (nrow < 0 || nrow >= Height) continue;
                for (int j = 0; j < MASK_WIDTH; ++j) {
                    int ncol = col - MASK_RADIUS + j;
                    if (ncol < 0 || ncol >= Width) continue;
                    sum += N[nrow * Width + ncol] *
                           h_Mck[i * MASK_WIDTH + j];
                }
            }
            P[row * Width + col] = sum;
        }
    }
}

/* 主机端 mask 副本 (已在 CPU 参考实现之前声明, 见上) */

int main(int argc, char **argv)
{
    int Width  = (argc > 1) ? atoi(argv[1]) : 2048;
    int Height = (argc > 2) ? atoi(argv[2]) : 2048;
    size_t nElem = (size_t)Width * Height;

    printf("=== 朴素 2D 卷积 (常量内存 mask) ===\n");
    printf("输入尺寸 %d x %d, mask %d x %d\n", Width, Height, MASK_WIDTH, MASK_WIDTH);

    /* 1. 主机端分配与初始化 */
    float *h_N = (float *)malloc(nElem * sizeof(float));
    float *h_P = (float *)malloc(nElem * sizeof(float));
    float *h_Pref = (float *)malloc(nElem * sizeof(float));
    if (!h_N || !h_P || !h_Pref) { fprintf(stderr, "host malloc failed\n"); return 1; }

    srand(1234);
    for (size_t i = 0; i < nElem; ++i) h_N[i] = (float)(rand() % 16);

    /* 2. 初始化 mask: 用经典的 5x5 高斯近似模糊核 */
    const float mk[MASK_WIDTH * MASK_WIDTH] = {
        1,  4,  6,  4, 1,
        4, 16, 24, 16, 4,
        6, 24, 36, 24, 6,
        4, 16, 24, 16, 4,
        1,  4,  6,  4, 1
    };
    float msum = 0.0f;
    for (int i = 0; i < MASK_WIDTH * MASK_WIDTH; ++i) msum += mk[i];
    for (int i = 0; i < MASK_WIDTH * MASK_WIDTH; ++i) {
        h_Mck[i] = mk[i] / msum;              /* 归一化, 保证模糊不改变整体亮度 */
    }

    /* 3. 把 mask 拷进常量内存: 必须用 cudaMemcpyToSymbol, 不能用 cudaMemcpy */
    CUDA_CHECK(cudaMemcpyToSymbol(Mc, h_Mck,
                                  MASK_WIDTH * MASK_WIDTH * sizeof(float)));

    /* 4. 设备端分配 + 传输输入 */
    float *d_N = NULL, *d_P = NULL;
    CUDA_CHECK(cudaMalloc((void **)&d_N, nElem * sizeof(float)));
    CUDA_CHECK(cudaMalloc((void **)&d_P, nElem * sizeof(float)));
    CUDA_CHECK(cudaMemcpy(d_N, h_N, nElem * sizeof(float), cudaMemcpyHostToDevice));

    /* 5. 启动配置: 16x16 的二维 block */
    dim3 block(16, 16);
    dim3 grid((Width + block.x - 1) / block.x, (Height + block.y - 1) / block.y);
    printf("grid = (%u, %u), block = (%u, %u) = %u 线程/block, %u warp/block\n",
           grid.x, grid.y, block.x, block.y, block.x * block.y,
           (block.x * block.y + 31) / 32);

    /* 6. 预热 + 计时 (5 次取平均, 用 cudaEvent 计时) */
    conv2d_naive_kernel<<<grid, block>>>(d_N, d_P, Width, Height);
    CUDA_CHECK(cudaGetLastError());
    CUDA_CHECK(cudaDeviceSynchronize());

    cudaEvent_t t0, t1;
    CUDA_CHECK(cudaEventCreate(&t0));
    CUDA_CHECK(cudaEventCreate(&t1));
    const int REPEAT = 5;
    CUDA_CHECK(cudaEventRecord(t0));
    for (int r = 0; r < REPEAT; ++r)
        conv2d_naive_kernel<<<grid, block>>>(d_N, d_P, Width, Height);
    CUDA_CHECK(cudaEventRecord(t1));
    CUDA_CHECK(cudaEventSynchronize(t1));
    CUDA_CHECK(cudaGetLastError());

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

    /* 7. 性能指标: FLOP 数 = 2 * MASK_WIDTH^2 * 输出元素数 */
    double flops  = 2.0 * MASK_WIDTH * MASK_WIDTH * (double)nElem;
    double bytes  = 2.0 * (double)nElem * sizeof(float);   /* 至少读写一遍输入/输出 */
    printf("耗时        : %.3f ms\n", ms);
    printf("GFLOP/s     : %.2f\n", flops / (ms * 1e6));
    printf("有效带宽     : %.2f GB/s (按 (读入+写出) 计算)\n",
           bytes / (ms * 1e6));

    /* 8. 取回结果并与 CPU 参考实现对比 */
    CUDA_CHECK(cudaMemcpy(h_P, d_P, nElem * sizeof(float), cudaMemcpyDeviceToHost));
    conv2d_cpu_reference(h_N, h_Pref, Width, Height);

    double maxErr = 0.0;
    for (size_t i = 0; i < nElem; ++i) {
        double e = fabs((double)h_P[i] - (double)h_Pref[i]);
        if (e > maxErr) maxErr = e;
    }
    printf("最大绝对误差 : %.3e\n", maxErr);
    printf("验证结果     : %s\n", (maxErr < 1e-3) ? "PASS" : "FAIL");

    CUDA_CHECK(cudaEventDestroy(t0));
    CUDA_CHECK(cudaEventDestroy(t1));
    CUDA_CHECK(cudaFree(d_N));
    CUDA_CHECK(cudaFree(d_P));
    free(h_N); free(h_P); free(h_Pref);
    return 0;
}

// 文件: conv2d_tiled_halo.cu
// 编译: nvcc -O3 -arch=sm_80 conv2d_tiled_halo.cu -o conv2d_tiled
// 运行: ./conv2d_tiled             # 默认 2048x2048 输入, 5x5 mask
//       ./conv2d_tiled 1024 1024
//
// 功能: 在同一份程序里实现并对比三个版本:
//       (A) conv2d_naive_kernel        朴素版 (对照, 与示例 1 相同)
//       (B) conv2d_tiled_kernel        分块版: block=(20x20), 共享内存 20x20 含 halo,
//                                      每线程搬 1 个元素, 前 256 个线程各算 1 个输出
//       (C) conv2d_tiled_coarse_kernel 分块+线程粗化: block=(8x8), 每线程算 2x2=4 个输出
//       三者输出全部与 CPU 参考实现比对, 并打印加速比

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

#define TILE_WIDTH   16
#define MASK_WIDTH   5
#define MASK_RADIUS  (MASK_WIDTH / 2)
#define IN_TILE      (TILE_WIDTH + MASK_WIDTH - 1)   /* 20: 输入 tile 的一边 */
#define COARSE       2                               /* 粗化: 每线程 2x2 个输出 */

static __constant__ float Mc[MASK_WIDTH * MASK_WIDTH];
static float h_Mck[MASK_WIDTH * MASK_WIDTH];

#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)

/* ---------------- (A) 朴素版: 一个线程一个输出, 边界分支 ---------------- */
__global__ void conv2d_naive_kernel(const float * __restrict__ N,
                                    float * __restrict__ P,
                                    int Width, int Height)
{
    int col = blockIdx.x * blockDim.x + threadIdx.x;
    int row = blockIdx.y * blockDim.y + threadIdx.y;
    if (row >= Height || col >= Width) return;

    float Pvalue = 0.0f;
    #pragma unroll
    for (int i = 0; i < MASK_WIDTH; ++i) {
        int nrow = row - MASK_RADIUS + i;
        if (nrow < 0 || nrow >= Height) continue;
        #pragma unroll
        for (int j = 0; j < MASK_WIDTH; ++j) {
            int ncol = col - MASK_RADIUS + j;
            if (ncol < 0 || ncol >= Width) continue;
            Pvalue += N[nrow * Width + ncol] * Mc[i * MASK_WIDTH + j];
        }
    }
    P[row * Width + col] = Pvalue;
}

/* ---------------- (B) 分块版: halo 协作加载 + 无分支计算 ---------------- */
__global__ void conv2d_tiled_kernel(const float * __restrict__ N,
                                    float * __restrict__ P,
                                    int Width, int Height)
{
    __shared__ float Ns[IN_TILE][IN_TILE];        /* 20 x 20 = 400 float = 1600 B */

    const int tx = threadIdx.x;                   /* 0..19 */
    const int ty = threadIdx.y;                   /* 0..19 */

    /* 该线程负责的输入坐标 (block 的左上角是输出瓦片左上角再向左上各推 radius) */
    const int row_i = blockIdx.y * TILE_WIDTH - MASK_RADIUS + ty;
    const int col_i = blockIdx.x * TILE_WIDTH - MASK_RADIUS + tx;

    /* 1. 协作加载 (TILE+2r)^2 个元素, 越界元素一律填 0 (zero padding) */
    float v = 0.0f;
    if (row_i >= 0 && row_i < Height && col_i >= 0 && col_i < Width)
        v = N[row_i * Width + col_i];
    Ns[ty][tx] = v;                               /* 分支之外: 所有线程都参与 */

    __syncthreads();                              /* 等整块 tile 就绪, 必须在分支外 */

    /* 2. 前 TILE x TILE 个线程各算一个输出 (计算阶段无任何分支) */
    if (ty < TILE_WIDTH && tx < TILE_WIDTH) {
        const int row_o = blockIdx.y * TILE_WIDTH + ty;
        const int col_o = blockIdx.x * TILE_WIDTH + tx;
        if (row_o < Height && col_o < Width) {
            float Pvalue = 0.0f;
            #pragma unroll
            for (int i = 0; i < MASK_WIDTH; ++i)
                #pragma unroll
                for (int j = 0; j < MASK_WIDTH; ++j)
                    Pvalue += Ns[ty + i][tx + j] * Mc[i * MASK_WIDTH + j];
            P[row_o * Width + col_o] = Pvalue;
        }
    }
}

/* ---------------- (C) 分块 + 线程粗化: 每线程 2x2 = 4 个输出 ---------------- */
__global__ void conv2d_tiled_coarse_kernel(const float * __restrict__ N,
                                           float * __restrict__ P,
                                           int Width, int Height)
{
    __shared__ float Ns[IN_TILE][IN_TILE];

    const int tx  = threadIdx.x;                  /* 0..7 */
    const int ty  = threadIdx.y;                  /* 0..7 */
    const int tid = ty * blockDim.x + tx;         /* 0..63 */
    const int nth = blockDim.x * blockDim.y;      /* 64 */

    const int row_b = blockIdx.y * TILE_WIDTH - MASK_RADIUS;
    const int col_b = blockIdx.x * TILE_WIDTH - MASK_RADIUS;

    /* 1. 64 个线程协作搬 400 个元素 (每个线程 6-7 次, 循环步长 = 线程数) */
    for (int l = tid; l < IN_TILE * IN_TILE; l += nth) {
        const int ly = l / IN_TILE;
        const int lx = l - ly * IN_TILE;
        const int gr = row_b + ly;
        const int gc = col_b + lx;
        float v = 0.0f;
        if (gr >= 0 && gr < Height && gc >= 0 && gc < Width)
            v = N[gr * Width + gc];
        Ns[ly][lx] = v;
    }
    __syncthreads();

    /* 2. 每个线程算 2x2 个输出, 4 个累加器常驻寄存器 */
    float acc[COARSE][COARSE];
    #pragma unroll
    for (int dy = 0; dy < COARSE; ++dy)
        #pragma unroll
        for (int dx = 0; dx < COARSE; ++dx)
            acc[dy][dx] = 0.0f;

    const int oy0 = ty * COARSE;                  /* tile 内的局部输出行 */
    const int ox0 = tx * COARSE;                  /* tile 内的局部输出列 */

    #pragma unroll
    for (int i = 0; i < MASK_WIDTH; ++i) {
        #pragma unroll
        for (int j = 0; j < MASK_WIDTH; ++j) {
            const float m = Mc[i * MASK_WIDTH + j];   /* 1 次常量读喂 4 个 FMA */
            #pragma unroll
            for (int dy = 0; dy < COARSE; ++dy)
                #pragma unroll
                for (int dx = 0; dx < COARSE; ++dx)
                    acc[dy][dx] += Ns[oy0 + dy + i][ox0 + dx + j] * m;
        }
    }

    #pragma unroll
    for (int dy = 0; dy < COARSE; ++dy) {
        #pragma unroll
        for (int dx = 0; dx < COARSE; ++dx) {
            const int row_o = blockIdx.y * TILE_WIDTH + oy0 + dy;
            const int col_o = blockIdx.x * TILE_WIDTH + ox0 + dx;
            if (row_o < Height && col_o < Width)
                P[row_o * Width + col_o] = acc[dy][dx];
        }
    }
}

/* ---------------- CPU 参考实现 ---------------- */
static void conv2d_cpu_reference(const float *N, float *P, int Width, int Height)
{
    for (int row = 0; row < Height; ++row) {
        for (int col = 0; col < Width; ++col) {
            float sum = 0.0f;
            for (int i = 0; i < MASK_WIDTH; ++i) {
                int nrow = row - MASK_RADIUS + i;
                if (nrow < 0 || nrow >= Height) continue;
                for (int j = 0; j < MASK_WIDTH; ++j) {
                    int ncol = col - MASK_RADIUS + j;
                    if (ncol < 0 || ncol >= Width) continue;
                    sum += N[nrow * Width + ncol] * h_Mck[i * MASK_WIDTH + j];
                }
            }
            P[row * Width + col] = sum;
        }
    }
}

/* 计时辅助: 通过函数指针启动 kernel 必须用 cudaLaunchKernel
   (<<<>>> 语法不能直接作用于函数指针) */
static float time_kernel(void (*kern)(const float *, float *, int, int),
                         dim3 grid, dim3 block, const float *d_N, float *d_P,
                         int W, int H, int repeat)
{
    void *args[] = { (void *)&d_N, (void *)&d_P, (void *)&W, (void *)&H };
    CUDA_CHECK(cudaLaunchKernel((const void *)kern, grid, block, args, 0, 0));
    CUDA_CHECK(cudaDeviceSynchronize());

    cudaEvent_t t0, t1;
    CUDA_CHECK(cudaEventCreate(&t0));
    CUDA_CHECK(cudaEventCreate(&t1));
    CUDA_CHECK(cudaEventRecord(t0));
    for (int r = 0; r < repeat; ++r)
        CUDA_CHECK(cudaLaunchKernel((const void *)kern, grid, block, args, 0, 0));
    CUDA_CHECK(cudaEventRecord(t1));
    CUDA_CHECK(cudaEventSynchronize(t1));
    CUDA_CHECK(cudaGetLastError());
    float ms = 0.0f;
    CUDA_CHECK(cudaEventElapsedTime(&ms, t0, t1));
    CUDA_CHECK(cudaEventDestroy(t0));
    CUDA_CHECK(cudaEventDestroy(t1));
    return ms / repeat;
}

static int verify(const float *h_P, const float *h_Pref, size_t n, const char *tag)
{
    double maxErr = 0.0;
    for (size_t i = 0; i < n; ++i) {
        double e = fabs((double)h_P[i] - (double)h_Pref[i]);
        if (e > maxErr) maxErr = e;
    }
    printf("  [%s] 最大绝对误差 %.3e -> %s\n", tag, maxErr,
           (maxErr < 1e-3) ? "PASS" : "FAIL");
    return (maxErr < 1e-3) ? 1 : 0;
}

int main(int argc, char **argv)
{
    int Width  = (argc > 1) ? atoi(argv[1]) : 2048;
    int Height = (argc > 2) ? atoi(argv[2]) : 2048;
    size_t nElem = (size_t)Width * Height;

    printf("=== 分块 2D 卷积 (halo 加载 + 线程粗化) ===\n");
    printf("输入 %d x %d, mask %d x %d, TILE %d, 输入 tile %d\n",
           Width, Height, MASK_WIDTH, MASK_WIDTH, TILE_WIDTH, IN_TILE);

    float *h_N    = (float *)malloc(nElem * sizeof(float));
    float *h_P    = (float *)malloc(nElem * sizeof(float));
    float *h_Pref = (float *)malloc(nElem * sizeof(float));
    if (!h_N || !h_P || !h_Pref) { fprintf(stderr, "host malloc failed\n"); return 1; }

    srand(1234);
    for (size_t i = 0; i < nElem; ++i) h_N[i] = (float)(rand() % 16);

    const float mk[MASK_WIDTH * MASK_WIDTH] = {
        1,  4,  6,  4, 1,
        4, 16, 24, 16, 4,
        6, 24, 36, 24, 6,
        4, 16, 24, 16, 4,
        1,  4,  6,  4, 1
    };
    float msum = 0.0f;
    for (int i = 0; i < MASK_WIDTH * MASK_WIDTH; ++i) msum += mk[i];
    for (int i = 0; i < MASK_WIDTH * MASK_WIDTH; ++i) h_Mck[i] = mk[i] / msum;

    CUDA_CHECK(cudaMemcpyToSymbol(Mc, h_Mck,
                                  MASK_WIDTH * MASK_WIDTH * sizeof(float)));

    float *d_N = NULL, *d_P = NULL;
    CUDA_CHECK(cudaMalloc((void **)&d_N, nElem * sizeof(float)));
    CUDA_CHECK(cudaMalloc((void **)&d_P, nElem * sizeof(float)));
    CUDA_CHECK(cudaMemcpy(d_N, h_N, nElem * sizeof(float), cudaMemcpyHostToDevice));

    const int REPEAT = 5;
    const double flops = 2.0 * MASK_WIDTH * MASK_WIDTH * (double)nElem;
    const double bytes = 2.0 * (double)nElem * sizeof(float);

    /* --- (A) 朴素版: block 16x16 --- */
    dim3 blkA(16, 16);
    dim3 grdA((Width + 15) / 16, (Height + 15) / 16);
    float msA = time_kernel(conv2d_naive_kernel, grdA, blkA, d_N, d_P, Width, Height, REPEAT);
    CUDA_CHECK(cudaMemcpy(h_P, d_P, nElem * sizeof(float), cudaMemcpyDeviceToHost));
    conv2d_cpu_reference(h_N, h_Pref, Width, Height);
    printf("\n(A) 朴素版        block(16,16)=256 线程 grid(%u,%u)\n", grdA.x, grdA.y);
    printf("    耗时 %.3f ms   GFLOP/s %.2f   有效带宽 %.2f GB/s\n",
           msA, flops / (msA * 1e6), bytes / (msA * 1e6));
    verify(h_P, h_Pref, nElem, "A");

    /* --- (B) 分块版: block (TILE+2r)x(TILE+2r) = 20x20 = 400 线程 --- */
    dim3 blkB(IN_TILE, IN_TILE);
    dim3 grdB((Width + TILE_WIDTH - 1) / TILE_WIDTH,
              (Height + TILE_WIDTH - 1) / TILE_WIDTH);
    float msB = time_kernel(conv2d_tiled_kernel, grdB, blkB, d_N, d_P, Width, Height, REPEAT);
    CUDA_CHECK(cudaMemcpy(h_P, d_P, nElem * sizeof(float), cudaMemcpyDeviceToHost));
    printf("\n(B) 分块版        block(%d,%d)=%d 线程 grid(%u,%u) 共享内存 %zu B/block\n",
           blkB.x, blkB.y, blkB.x * blkB.y, grdB.x, grdB.y,
           sizeof(float) * IN_TILE * IN_TILE);
    printf("    耗时 %.3f ms   GFLOP/s %.2f   有效带宽 %.2f GB/s   相对 A 加速 %.2fx\n",
           msB, flops / (msB * 1e6), bytes / (msB * 1e6), msA / msB);
    verify(h_P, h_Pref, nElem, "B");

    /* --- (C) 分块 + 粗化版: block 8x8 = 64 线程, 每线程 2x2 输出 --- */
    dim3 blkC(TILE_WIDTH / COARSE, TILE_WIDTH / COARSE);
    dim3 grdC((Width + TILE_WIDTH - 1) / TILE_WIDTH,
              (Height + TILE_WIDTH - 1) / TILE_WIDTH);
    float msC = time_kernel(conv2d_tiled_coarse_kernel, grdC, blkC, d_N, d_P, Width, Height, REPEAT);
    CUDA_CHECK(cudaMemcpy(h_P, d_P, nElem * sizeof(float), cudaMemcpyDeviceToHost));
    printf("\n(C) 分块+粗化版    block(%d,%d)=%d 线程, 每线程 %dx%d=4 个输出\n",
           blkC.x, blkC.y, blkC.x * blkC.y, COARSE, COARSE);
    printf("    耗时 %.3f ms   GFLOP/s %.2f   有效带宽 %.2f GB/s   相对 A 加速 %.2fx, 相对 B %.2fx\n",
           msC, flops / (msC * 1e6), bytes / (msC * 1e6), msA / msC, msB / msC);
    verify(h_P, h_Pref, nElem, "C");

    printf("\n理论分析: 分块把每输出全局读从 25 次降到 (TILE+2r)^2/TILE^2 = %.4f 次\n",
           (double)(IN_TILE * IN_TILE) / (double)(TILE_WIDTH * TILE_WIDTH));
    printf("          全局内存访问减少倍数 = TILE^2*MW^2/(TILE+MW-1)^2 = %.2fx\n",
           (double)(TILE_WIDTH * TILE_WIDTH * MASK_WIDTH * MASK_WIDTH) /
           (double)(IN_TILE * IN_TILE));

    CUDA_CHECK(cudaFree(d_N));
    CUDA_CHECK(cudaFree(d_P));
    free(h_N); free(h_P); free(h_Pref);
    return 0;
}

// 文件: conv3d_tiled_coarse.cu
// 编译: nvcc -O3 -arch=sm_80 conv3d_tiled_coarse.cu -o conv3d_tiled
// 运行: ./conv3d_tiled            # 默认 64x64x64 输入, 5x5x5 mask
//       ./conv3d_tiled 96 96 96
//
// 功能: 3D 卷积的三个版本对比 (对应 Lab 4: 3D Convolution / MP-4)
//       (A) conv3d_naive_kernel          朴素版: 一个线程一个输出, 每输出 125 次全局读
//       (B) conv3d_tiled_kernel          分块版: 共享内存 12x12x12 含 halo,
//                                        每个线程计算 COARSE_Z = 2 个 z 方向相邻输出
//       (C) conv3d_pitched_naive_kernel  用 cudaMalloc3D / cudaMemcpy3D /
//                                        cudaPitchedPtr 组织数据, 演示 pitch 寻址
//       kernel 用一维展平索引 idx = x + y*W + z*W*H (三维数组的线性化);
//       支持两种输出约定:
//         origin = 0 -> centered 形式, 输出尺寸 = 输入尺寸, 边界补零 (zero padding)
//         origin = r -> Lab 约定, P[z][y][x] = Σ M[i][j][k]*N[z+i][y+j][x+k],
//                       输出尺寸 = (W-M+1) x (H-M+1) x (D-M+1)

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

#define MASK_W   5
#define R        (MASK_W / 2)          /* 2 */

#define TX       8                     /* block 的 x 维线程数 */
#define TY       8                     /* block 的 y 维线程数 */
#define TZ       4                     /* block 的 z 维线程数 (仅 z 方向粗化) */
#define COARSE_Z 2                     /* 每个线程计算 COARSE_Z 个相邻 z 输出 */
#define TILE_X   TX
#define TILE_Y   TY
#define TILE_Z   (TZ * COARSE_Z)       /* block 负责的 z 方向输出数 = 8 */
#define SX       (TILE_X + MASK_W - 1) /* 12 */
#define SY       (TILE_Y + MASK_W - 1) /* 12 */
#define SZ       (TILE_Z + MASK_W - 1) /* 12 */
#define PAD      0                     /* 改成 1 则把 Ns 的 x 行长补到 13, 见 bank 分析 */

static __constant__ float Mc3[MASK_W][MASK_W][MASK_W];   /* 125*4 = 500 B << 64 KB */

#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)

/* ---------------- (A) 朴素 3D 卷积: 一个线程一个输出 ---------------- */
__global__ void conv3d_naive_kernel(const float * __restrict__ N,
                                    float * __restrict__ P,
                                    int W, int H, int D,
                                    int Ws, int Hs, int Ds, int origin)
{
    const int ox = blockIdx.x * blockDim.x + threadIdx.x;
    const int oy = blockIdx.y * blockDim.y + threadIdx.y;
    const int oz = blockIdx.z * blockDim.z + threadIdx.z;
    if (ox >= Ws \|\| oy >= Hs \|\| oz >= Ds) return;

    float sum = 0.0f;
    #pragma unroll
    for (int i = 0; i < MASK_W; ++i) {
        const int gz = oz + origin - R + i;
        if (gz < 0 \|\| gz >= D) continue;
        #pragma unroll
        for (int j = 0; j < MASK_W; ++j) {
            const int gy = oy + origin - R + j;
            if (gy < 0 \|\| gy >= H) continue;
            #pragma unroll
            for (int k = 0; k < MASK_W; ++k) {
                const int gx = ox + origin - R + k;
                if (gx < 0 \|\| gx >= W) continue;
                sum += N[(gz * H + gy) * W + gx] * Mc3[i][j][k];
            }
        }
    }
    P[(oz * Hs + oy) * Ws + ox] = sum;
}

/* ---------------- (B) 分块 3D 卷积: halo 加载 + z 方向粗化 ---------------- */
__global__ void conv3d_tiled_kernel(const float * __restrict__ N,
                                    float * __restrict__ P,
                                    int W, int H, int D,
                                    int Ws, int Hs, int Ds, int origin)
{
#if PAD
    __shared__ float Ns[SZ][SY][SX + 1];          /* 补 1 列, 行长 13 (奇数) */
#else
    __shared__ float Ns[SZ][SY][SX];              /* 12 x 12 x 12 = 6912 B */
#endif

    const int tx = threadIdx.x;                   /* 0..7  */
    const int ty = threadIdx.y;                   /* 0..7  */
    const int tz = threadIdx.z;                   /* 0..3  */
    const int tid = (tz * blockDim.y + ty) * blockDim.x + tx;
    const int nth = blockDim.x * blockDim.y * blockDim.z;   /* 256 */

    /* block 在输入空间中的起点 (origin 决定 centered / Lab 两种约定) */
    const int bx = blockIdx.x * TILE_X + origin - R;
    const int by = blockIdx.y * TILE_Y + origin - R;
    const int bz = blockIdx.z * TILE_Z + origin - R;

    /* 1. 256 个线程协作搬 12*12*12 = 1728 个元素进共享内存 (每线程 6-7 次)
          越界元素补 0 —— 这就是 halo 的 zero padding */
    for (int l = tid; l < SX * SY * SZ; l += nth) {
        const int lx = l % SX;
        const int ly = (l / SX) % SY;
        const int lz = l / (SX * SY);
        const int gx = bx + lx, gy = by + ly, gz = bz + lz;
        float v = 0.0f;
        if (gx >= 0 && gx < W && gy >= 0 && gy < H && gz >= 0 && gz < D)
            v = N[(gz * H + gy) * W + gx];        /* 一维展平索引 */
        Ns[lz][ly][lx] = v;
    }
    __syncthreads();

    /* 2. 每个线程算 COARSE_Z 个 z 方向相邻的输出 */
    float acc[COARSE_Z];
    #pragma unroll
    for (int c = 0; c < COARSE_Z; ++c) acc[c] = 0.0f;

    #pragma unroll
    for (int i = 0; i < MASK_W; ++i) {
        #pragma unroll
        for (int j = 0; j < MASK_W; ++j) {
            #pragma unroll
            for (int k = 0; k < MASK_W; ++k) {
                const float m = Mc3[i][j][k];     /* 1 次常量读 */
                #pragma unroll
                for (int c = 0; c < COARSE_Z; ++c)  /* 喂给 2 个 FMA */
                    acc[c] += Ns[tz * COARSE_Z + c + i][ty + j][tx + k] * m;
            }
        }
    }

    /* 3. 写回 COARSE_Z 个输出 (三维线性化: idx = x + y*Ws + z*Ws*Hs) */
    #pragma unroll
    for (int c = 0; c < COARSE_Z; ++c) {
        const int ox = blockIdx.x * TILE_X + tx;
        const int oy = blockIdx.y * TILE_Y + ty;
        const int oz = blockIdx.z * TILE_Z + tz * COARSE_Z + c;
        if (ox < Ws && oy < Hs && oz < Ds)
            P[(oz * Hs + oy) * Ws + ox] = acc[c];
    }
}

/* ---------------- (C) 用 cudaPitchedPtr 的朴素 kernel (演示 pitch 索引) -------
   陷阱提醒: 输入与输出是两次独立的 cudaMalloc3D, 它们的 pitch 可能不同,
   因此必须分别传入 pitchN 与 pitchP, 换算成 float 行步长 Wpn / Wpp 使用。      */
__global__ void conv3d_pitched_naive_kernel(const float * __restrict__ Np,
                                            float * __restrict__ Pp,
                                            size_t pitchNBytes,
                                            size_t pitchPBytes,
                                            int W, int H, int D,
                                            int Ws, int Hs, int Ds, int origin)
{
    const int Wpn = (int)(pitchNBytes / sizeof(float));  /* 输入行内 float 数 >= W  */
    const int Wpp = (int)(pitchPBytes / sizeof(float));  /* 输出行内 float 数 >= Ws */
    const int ox = blockIdx.x * blockDim.x + threadIdx.x;
    const int oy = blockIdx.y * blockDim.y + threadIdx.y;
    const int oz = blockIdx.z * blockDim.z + threadIdx.z;
    if (ox >= Ws \|\| oy >= Hs \|\| oz >= Ds) return;

    float sum = 0.0f;
    #pragma unroll
    for (int i = 0; i < MASK_W; ++i) {
        const int gz = oz + origin - R + i;
        if (gz < 0 \|\| gz >= D) continue;
        #pragma unroll
        for (int j = 0; j < MASK_W; ++j) {
            const int gy = oy + origin - R + j;
            if (gy < 0 \|\| gy >= H) continue;
            #pragma unroll
            for (int k = 0; k < MASK_W; ++k) {
                const int gx = ox + origin - R + k;
                if (gx < 0 \|\| gx >= W) continue;
                sum += Np[gz * Wpn * H + gy * Wpn + gx] * Mc3[i][j][k];
            }
        }
    }
    Pp[oz * Wpp * Hs + oy * Wpp + ox] = sum;     /* 输出坐标 oz/oy/ox, 行步长 Wpp */
}

/* ---------------- CPU 参考实现 (同时支持两种约定) ---------------- */
static float h_Mc3[MASK_W][MASK_W][MASK_W];

static void conv3d_cpu_reference(const float *N, float *P,
                                 int W, int H, int D,
                                 int Ws, int Hs, int Ds, int origin)
{
    for (int oz = 0; oz < Ds; ++oz)
        for (int oy = 0; oy < Hs; ++oy)
            for (int ox = 0; ox < Ws; ++ox) {
                float sum = 0.0f;
                for (int i = 0; i < MASK_W; ++i) {
                    int gz = oz + origin - R + i;
                    if (gz < 0 \|\| gz >= D) continue;
                    for (int j = 0; j < MASK_W; ++j) {
                        int gy = oy + origin - R + j;
                        if (gy < 0 \|\| gy >= H) continue;
                        for (int k = 0; k < MASK_W; ++k) {
                            int gx = ox + origin - R + k;
                            if (gx < 0 \|\| gx >= W) continue;
                            sum += N[(gz * H + gy) * W + gx] * h_Mc3[i][j][k];
                        }
                    }
                }
                P[(oz * Hs + oy) * Ws + ox] = sum;
            }
}

/* h_Mc3 (主机端 mask 副本) 已在 CPU 参考实现之前声明 */

static float time_flat_kernel(void (*kern)(const float *, float *, int, int, int,
                                           int, int, int, int),
                              dim3 grid, dim3 block,
                              const float *d_N, float *d_P,
                              int W, int H, int D, int Ws, int Hs, int Ds,
                              int origin, int repeat)
{
    void *args[] = { (void *)&d_N, (void *)&d_P, (void *)&W, (void *)&H,
                     (void *)&D, (void *)&Ws, (void *)&Hs, (void *)&Ds,
                     (void *)&origin };
    CUDA_CHECK(cudaLaunchKernel((const void *)kern, grid, block, args, 0, 0));
    CUDA_CHECK(cudaDeviceSynchronize());
    cudaEvent_t t0, t1;
    CUDA_CHECK(cudaEventCreate(&t0));
    CUDA_CHECK(cudaEventCreate(&t1));
    CUDA_CHECK(cudaEventRecord(t0));
    for (int r = 0; r < repeat; ++r)
        CUDA_CHECK(cudaLaunchKernel((const void *)kern, grid, block, args, 0, 0));
    CUDA_CHECK(cudaEventRecord(t1));
    CUDA_CHECK(cudaEventSynchronize(t1));
    CUDA_CHECK(cudaGetLastError());
    float ms = 0.0f;
    CUDA_CHECK(cudaEventElapsedTime(&ms, t0, t1));
    CUDA_CHECK(cudaEventDestroy(t0));
    CUDA_CHECK(cudaEventDestroy(t1));
    return ms / repeat;
}

static void verify3d(const float *a, const float *b, size_t n, const char *tag)
{
    double maxErr = 0.0;
    for (size_t i = 0; i < n; ++i) {
        double e = fabs((double)a[i] - (double)b[i]);
        if (e > maxErr) maxErr = e;
    }
    printf("  [%s] 最大绝对误差 %.3e -> %s\n", tag, maxErr,
           (maxErr < 1e-3) ? "PASS" : "FAIL");
}

int main(int argc, char **argv)
{
    const int W = (argc > 1) ? atoi(argv[1]) : 64;
    const int H = (argc > 2) ? atoi(argv[2]) : 64;
    const int D = (argc > 3) ? atoi(argv[3]) : 64;
    const size_t nIn = (size_t)W * H * D;

    printf("=== 3D 分块卷积 (Lab 4: 3D Convolution / MP-4) ===\n");
    printf("输入 %d x %d x %d, mask %dx%dx%d, 共享内存 tile %dx%dx%d\n",
           W, H, D, MASK_W, MASK_W, MASK_W, SZ, SY, SX);
    printf("block = (%d,%d,%d) = %d 线程, 每线程 %d 个 z 输出, 每 block 输出 %d 个\n",
           TX, TY, TZ, TX * TY * TZ, COARSE_Z, TILE_X * TILE_Y * TILE_Z);

    float *h_N = (float *)malloc(nIn * sizeof(float));
    if (!h_N) { fprintf(stderr, "host malloc failed\n"); return 1; }
    srand(4321);
    for (size_t i = 0; i < nIn; ++i) h_N[i] = (float)(rand() % 8);

    /* 3D mask: 用 5^3 的分离式高斯核 (可分离核能保证数值温和) */
    static const float g1[MASK_W] = {1, 4, 6, 4, 1};
    double gs = 0.0;
    for (int i = 0; i < MASK_W; ++i)
        for (int j = 0; j < MASK_W; ++j)
            for (int k = 0; k < MASK_W; ++k) {
                h_Mc3[i][j][k] = g1[i] * g1[j] * g1[k];
                gs += h_Mc3[i][j][k];
            }
    for (int i = 0; i < MASK_W; ++i)
        for (int j = 0; j < MASK_W; ++j)
            for (int k = 0; k < MASK_W; ++k)
                h_Mc3[i][j][k] = (float)(h_Mc3[i][j][k] / gs);

    CUDA_CHECK(cudaMemcpyToSymbol(Mc3, h_Mc3, sizeof(h_Mc3)));

    float *d_N = NULL, *d_P = NULL;
    CUDA_CHECK(cudaMalloc((void **)&d_N, nIn * sizeof(float)));
    CUDA_CHECK(cudaMalloc((void **)&d_P, nIn * sizeof(float)));   /* 最大输出尺寸 */
    CUDA_CHECK(cudaMemcpy(d_N, h_N, nIn * sizeof(float), cudaMemcpyHostToDevice));

    float *h_P    = (float *)malloc(nIn * sizeof(float));
    float *h_Pref = (float *)malloc(nIn * sizeof(float));
    if (!h_P \|\| !h_Pref) { fprintf(stderr, "host malloc failed\n"); return 1; }

    const int REPEAT = 5;

    /* ======================= 约定一: Lab 约定 (origin = R) ======================= */
    /* 输出 (W-M+1) x (H-M+1) x (D-M+1), P[z][y][x] = Σ M[i][j][k]*N[z+i][y+j][x+k] */
    {
        const int Ws = W - MASK_W + 1, Hs = H - MASK_W + 1, Ds = D - MASK_W + 1;
        const int origin = R;
        const size_t nOut = (size_t)Ws * Hs * Ds;
        const double flops = 2.0 * MASK_W * MASK_W * MASK_W * (double)nOut;
        const double bytes = ((double)nIn + (double)nOut) * sizeof(float);

        dim3 blkA(8, 8, 4);
        dim3 grdA((Ws + 7) / 8, (Hs + 7) / 8, (Ds + 3) / 4);
        float msA = time_flat_kernel(conv3d_naive_kernel, grdA, blkA, d_N, d_P,
                                     W, H, D, Ws, Hs, Ds, origin, REPEAT);
        CUDA_CHECK(cudaMemcpy(h_P, d_P, nOut * sizeof(float), cudaMemcpyDeviceToHost));
        conv3d_cpu_reference(h_N, h_Pref, W, H, D, Ws, Hs, Ds, origin);
        printf("\n[Lab 约定] 输出 %d x %d x %d\n", Ws, Hs, Ds);
        printf("  (A) 朴素 3D     %.3f ms  %.2f GFLOP/s  %.2f GB/s\n",
               msA, flops / (msA * 1e6), bytes / (msA * 1e6));
        verify3d(h_P, h_Pref, nOut, "A");

        dim3 blkB(TX, TY, TZ);
        dim3 grdB((Ws + TILE_X - 1) / TILE_X,
                  (Hs + TILE_Y - 1) / TILE_Y,
                  (Ds + TILE_Z - 1) / TILE_Z);
        float msB = time_flat_kernel(conv3d_tiled_kernel, grdB, blkB, d_N, d_P,
                                     W, H, D, Ws, Hs, Ds, origin, REPEAT);
        CUDA_CHECK(cudaMemcpy(h_P, d_P, nOut * sizeof(float), cudaMemcpyDeviceToHost));
        printf("  (B) 分块 3D     %.3f ms  %.2f GFLOP/s  %.2f GB/s  加速 %.2fx\n",
               msB, flops / (msB * 1e6), bytes / (msB * 1e6), msA / msB);
        verify3d(h_P, h_Pref, nOut, "B");

        printf("  理论: 减少倍数 R = %d^3*%d^3/%d^3 = %.2fx, halo 放大 = %.3f\n",
               TILE_X, MASK_W, SX, (double)(TILE_X * TILE_Y * TILE_Z * MASK_W * MASK_W * MASK_W) /
               (double)(SX * SY * SZ), (double)(SX * SY * SZ) / (double)(TILE_X * TILE_Y * TILE_Z));
    }

    /* =================== 约定二: centered + 边界补零 (origin = 0) =================== */
    /* 输出尺寸 = 输入尺寸, 图像外一律视作 0 (zero padding) */
    {
        const int Ws = W, Hs = H, Ds = D;
        const int origin = 0;
        const size_t nOut = nIn;
        const double flops = 2.0 * MASK_W * MASK_W * MASK_W * (double)nOut;
        const double bytes = ((double)nIn + (double)nOut) * sizeof(float);

        dim3 blkB(TX, TY, TZ);
        dim3 grdB((Ws + TILE_X - 1) / TILE_X,
                  (Hs + TILE_Y - 1) / TILE_Y,
                  (Ds + TILE_Z - 1) / TILE_Z);
        float msB = time_flat_kernel(conv3d_tiled_kernel, grdB, blkB, d_N, d_P,
                                     W, H, D, Ws, Hs, Ds, origin, REPEAT);
        CUDA_CHECK(cudaMemcpy(h_P, d_P, nOut * sizeof(float), cudaMemcpyDeviceToHost));
        conv3d_cpu_reference(h_N, h_Pref, W, H, D, Ws, Hs, Ds, origin);
        printf("\n[centered + 补零] 输出 %d x %d x %d (与输入同尺寸)\n", Ws, Hs, Ds);
        printf("  (B) 分块 3D     %.3f ms  %.2f GFLOP/s  %.2f GB/s\n",
               msB, flops / (msB * 1e6), bytes / (msB * 1e6));
        verify3d(h_P, h_Pref, nOut, "B-padded");
    }

    /* ============ 约定三: cudaMalloc3D + cudaMemcpy3D + cudaPitchedPtr ============ */
    {
        const int Ws = W - MASK_W + 1, Hs = H - MASK_W + 1, Ds = D - MASK_W + 1;
        const int origin = R;
        const size_t nOut = (size_t)Ws * Hs * Ds;

        cudaExtent extent = make_cudaExtent((size_t)W * sizeof(float), H, D);
        cudaPitchedPtr d_Np;
        CUDA_CHECK(cudaMalloc3D(&d_Np, extent));
        cudaMemcpy3DParms p = {0};
        p.srcPtr = make_cudaPitchedPtr((void *)h_N, (size_t)W * sizeof(float), W, H);
        p.dstPtr = d_Np;
        p.extent = extent;
        p.kind   = cudaMemcpyHostToDevice;
        CUDA_CHECK(cudaMemcpy3D(&p));

        cudaExtent oext = make_cudaExtent((size_t)Ws * sizeof(float), Hs, Ds);
        cudaPitchedPtr d_Pp;
        CUDA_CHECK(cudaMalloc3D(&d_Pp, oext));

        dim3 blk(8, 8, 4);
        dim3 grd((Ws + 7) / 8, (Hs + 7) / 8, (Ds + 3) / 4);
        conv3d_pitched_naive_kernel<<<grd, blk>>>((const float *)d_Np.ptr,
                                                  (float *)d_Pp.ptr,
                                                  d_Np.pitch, d_Pp.pitch,
                                                  W, H, D, Ws, Hs, Ds, origin);
        CUDA_CHECK(cudaGetLastError());
        CUDA_CHECK(cudaDeviceSynchronize());

        /* 把 pitched 结果按行拷回紧凑数组后验证 */
        float *h_Pp = (float *)malloc((size_t)Ws * Hs * Ds * sizeof(float));
        if (!h_Pp) { fprintf(stderr, "host malloc failed\n"); return 1; }
        cudaMemcpy3DParms q = {0};
        q.srcPtr = d_Pp;
        q.dstPtr = make_cudaPitchedPtr((void *)h_Pp, (size_t)Ws * sizeof(float), Ws, Hs);
        q.extent = oext;
        q.kind   = cudaMemcpyDeviceToHost;
        CUDA_CHECK(cudaMemcpy3D(&q));

        conv3d_cpu_reference(h_N, h_Pref, W, H, D, Ws, Hs, Ds, origin);
        printf("\n[cudaPitchedPtr 版本]\n");
        printf("  d_Np.pitch = %zu 字节 (W*4 = %zu 字节, 硬件对齐后可能更大)\n",
               d_Np.pitch, (size_t)W * sizeof(float));
        printf("  d_Pp.pitch = %zu 字节 (Ws*4 = %zu 字节)\n",
               d_Pp.pitch, (size_t)Ws * sizeof(float));
        printf("  kernel 中行内索引必须用 Wpn = pitchN/4 = %zu (输入),"
               " Wpp = pitchP/4 = %zu (输出),\n", d_Np.pitch / sizeof(float),
               d_Pp.pitch / sizeof(float));
        printf("  而不能用输入的 W = %d 或输出的 Ws = %d\n", W, Ws);
        verify3d(h_Pp, h_Pref, nOut, "C-pitched");

        CUDA_CHECK(cudaFree(d_Np.ptr));
        CUDA_CHECK(cudaFree(d_Pp.ptr));
        free(h_Pp);
    }

    CUDA_CHECK(cudaFree(d_N));
    CUDA_CHECK(cudaFree(d_P));
    free(h_N); free(h_P); free(h_Pref);
    return 0;
}