Lecture 4: 并行模式二 —— 基础矩阵乘法及其性能瓶颈 (对应 Lab 2: Simple Matrix Multiply)

// 文件: matmul_naive.cu
// 主题: ECE408 Lecture 4 / Lab 2 -- 朴素矩阵乘法（每线程一个输出元素，不使用共享内存）
// 编译: nvcc -O3 -arch=sm_80 matmul_naive.cu -o matmul_naive
// 运行: ./matmul_naive 1024        （参数为矩阵边长 Width，建议 512 / 1024 / 2048）

#include <cstdio>
#include <cstdlib>
#include <cmath>
#include <ctime>
#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 BLOCK_WIDTH   16    // 16 x 16 = 256 threads = 8 warps per block
#define WARMUP_ITERS   5
#define TIMED_ITERS   20

// ============================== device kernel ==============================
// 与讲义中的 Simple Matrix Multiplication Kernel 完全一致
__global__ void matmulNaiveKernel(const float* __restrict__ d_M,
                                  const float* __restrict__ d_N,
                                  float* __restrict__ d_P,
                                  int Width)
{
    // Calculate the row index of the d_P element and d_M
    int Row = blockIdx.y * blockDim.y + threadIdx.y;
    // Calculate the column index of d_P and d_N
    int Col = blockIdx.x * blockDim.x + threadIdx.x;

    if (Row < Width && Col < Width) {
        float Pvalue = 0.0f;
        // each thread computes one element of the block sub-matrix
        for (int k = 0; k < Width; ++k) {
            Pvalue += d_M[Row * Width + k] * d_N[k * Width + Col];
        }
        d_P[Row * Width + Col] = Pvalue;
    }
}

// ============================== host utilities ==============================
static void initMatrix(float* a, int n, unsigned int seed)
{
    srand(seed);
    for (int i = 0; i < n * n; ++i) {
        a[i] = (float)(rand() % 7) * 0.25f;   // 取值 0.00 ~ 1.50，便于结果校验
    }
}

// CPU 参考实现：i-k-j 顺序，最内层沿行连续访问，比 i-j-k 快得多
static void cpuMatmul(const float* A, const float* B, float* C, int n)
{
    for (int i = 0; i < n * n; ++i) C[i] = 0.0f;
    for (int i = 0; i < n; ++i) {
        for (int k = 0; k < n; ++k) {
            const float  a    = A[i * n + k];
            const float* Brow = B + k * n;
            float*       Crow = C + i * n;
            for (int j = 0; j < n; ++j) {
                Crow[j] += a * Brow[j];
            }
        }
    }
}

// 每个 SM 的 FP32 单元数（估算理论峰值算力用）
static int fp32CoresPerSM(int major, int minor)
{
    if (major == 7) return 64;                        // Volta (sm_70) / Turing (sm_75)
    if (major == 8) return (minor == 0) ? 64 : 128;   // A100 (sm_80)=64, GA10x (sm_86/89)=128
    if (major >= 9) return 128;                       // Hopper (sm_90) 及以后
    return 128;
}

// ================================== main ===================================
int main(int argc, char** argv)
{
    const int    N     = (argc > 1) ? atoi(argv[1]) : 1024;
    const size_t bytes = (size_t)N * (size_t)N * sizeof(float);

    // ---------- 1. 设备属性：理论峰值带宽与峰值算力 ----------
    CUDA_CHECK(cudaSetDevice(0));
    cudaDeviceProp prop;
    CUDA_CHECK(cudaGetDeviceProperties(&prop, 0));

    const double peakBW   = 2.0 * (double)prop.memoryClockRate * 1.0e3
                          * ((double)prop.memoryBusWidth / 8.0) / 1.0e9;   // GB/s
    const double peakFP32 = (double)prop.multiProcessorCount
                          * (double)fp32CoresPerSM(prop.major, prop.minor)
                          * 2.0 * (double)prop.clockRate * 1.0e3 / 1.0e12; // TFLOPS

    printf("== Device ==\n");
    printf("  name                 : %s\n", prop.name);
    printf("  compute capability   : sm_%d%d\n", prop.major, prop.minor);
    printf("  SMs                  : %d\n", prop.multiProcessorCount);
    printf("  core clock           : %.3f GHz\n", prop.clockRate * 1e-6);
    printf("  memory clock         : %.3f GHz, bus %d bit\n",
           prop.memoryClockRate * 1e-6, prop.memoryBusWidth);
    printf("  peak bandwidth       : %.1f GB/s\n", peakBW);
    printf("  peak FP32            : %.2f TFLOPS\n", peakFP32);
    printf("  machine balance      : %.2f FLOP/Byte\n",
           peakFP32 * 1e12 / (peakBW * 1e9));
    printf("  L2 cache             : %.1f MB\n", prop.l2CacheSize / 1048576.0);
    printf("  max threads per SM   : %d\n", prop.maxThreadsPerMultiProcessor);

    // ---------- 2. 主机端数据与 CPU 参考结果 ----------
    float* h_M   = (float*)malloc(bytes);
    float* h_N   = (float*)malloc(bytes);
    float* h_P   = (float*)malloc(bytes);
    float* h_Ref = (float*)malloc(bytes);
    if (!h_M || !h_N || !h_P || !h_Ref) {
        fprintf(stderr, "host malloc failed\n");
        exit(EXIT_FAILURE);
    }
    initMatrix(h_M, N, 1u);
    initMatrix(h_N, N, 2u);

    const clock_t cpuT0 = clock();
    cpuMatmul(h_M, h_N, h_Ref, N);
    const double cpuSec = (double)(clock() - cpuT0) / CLOCKS_PER_SEC;

    // ---------- 3. 设备端分配与拷贝 ----------
    float *d_M = NULL, *d_N = NULL, *d_P = NULL;
    CUDA_CHECK(cudaMalloc((void**)&d_M, bytes));
    CUDA_CHECK(cudaMalloc((void**)&d_N, bytes));
    CUDA_CHECK(cudaMalloc((void**)&d_P, bytes));
    CUDA_CHECK(cudaMemcpy(d_M, h_M, bytes, cudaMemcpyHostToDevice));
    CUDA_CHECK(cudaMemcpy(d_N, h_N, bytes, cudaMemcpyHostToDevice));

    // ---------- 4. 执行配置（讲义原文的写法：向上取整） ----------
    dim3 dimBlock(BLOCK_WIDTH, BLOCK_WIDTH, 1);
    dim3 dimGrid((N + BLOCK_WIDTH - 1) / BLOCK_WIDTH,
                 (N + BLOCK_WIDTH - 1) / BLOCK_WIDTH, 1);
    printf("\n== Launch configuration ==\n");
    printf("  dimBlock = (%u, %u, 1) = %u threads = %u warps\n",
           dimBlock.x, dimBlock.y, dimBlock.x * dimBlock.y,
           dimBlock.x * dimBlock.y / 32);
    printf("  dimGrid  = (%u, %u, 1) = %u blocks\n",
           dimGrid.x, dimGrid.y, dimGrid.x * dimGrid.y);
    printf("  total threads        : %llu\n",
           (unsigned long long)dimGrid.x * dimGrid.y * dimBlock.x * dimBlock.y);

    // ---------- 5. 资源占用（寄存器 / 占用率） ----------
    cudaFuncAttributes attr;
    CUDA_CHECK(cudaFuncGetAttributes(&attr, matmulNaiveKernel));
    int blocksPerSM = 0;
    CUDA_CHECK(cudaOccupancyMaxActiveBlocksPerMultiprocessor(
        &blocksPerSM, matmulNaiveKernel, BLOCK_WIDTH * BLOCK_WIDTH, 0));
    const int threadsPerSM = blocksPerSM * BLOCK_WIDTH * BLOCK_WIDTH;
    printf("\n== Occupancy ==\n");
    printf("  registers per thread : %d\n", attr.numRegs);
    printf("  static shared memory : %zu Byte\n", attr.sharedSizeBytes);
    printf("  blocks per SM        : %d\n", blocksPerSM);
    printf("  threads per SM       : %d  -> occupancy %.1f%%\n",
           threadsPerSM,
           100.0 * threadsPerSM / (double)prop.maxThreadsPerMultiProcessor);

    // ---------- 6. 计时：先热身，再取多次平均 ----------
    for (int i = 0; i < WARMUP_ITERS; ++i) {
        matmulNaiveKernel<<<dimGrid, dimBlock>>>(d_M, d_N, d_P, N);
    }
    CUDA_CHECK(cudaGetLastError());
    CUDA_CHECK(cudaDeviceSynchronize());

    cudaEvent_t evStart, evStop;
    CUDA_CHECK(cudaEventCreate(&evStart));
    CUDA_CHECK(cudaEventCreate(&evStop));
    CUDA_CHECK(cudaEventRecord(evStart));
    for (int i = 0; i < TIMED_ITERS; ++i) {
        matmulNaiveKernel<<<dimGrid, dimBlock>>>(d_M, d_N, d_P, N);
    }
    CUDA_CHECK(cudaEventRecord(evStop));
    CUDA_CHECK(cudaEventSynchronize(evStop));
    float totalMs = 0.0f;
    CUDA_CHECK(cudaEventElapsedTime(&totalMs, evStart, evStop));
    const double ms = (double)totalMs / TIMED_ITERS;

    // ---------- 7. 结果校验 ----------
    CUDA_CHECK(cudaMemcpy(h_P, d_P, bytes, cudaMemcpyDeviceToHost));
    double maxAbsErr = 0.0, maxRef = 0.0;
    for (int i = 0; i < N * N; ++i) {
        const double e = fabs((double)h_P[i] - (double)h_Ref[i]);
        if (e > maxAbsErr) maxAbsErr = e;
        if (fabs((double)h_Ref[i]) > maxRef) maxRef = fabs((double)h_Ref[i]);
    }
    const double relErr = (maxRef > 0.0) ? (maxAbsErr / maxRef) : maxAbsErr;

    // ---------- 8. 性能、算术强度与 Roofline 上限 ----------
    const double flops       = 2.0 * (double)N * N * N;
    const double bytesNeeded = 4.0 * flops;                     // 4 Byte per FLOP
    const double gflops      = flops / (ms * 1.0e-3) / 1.0e9;
    const double reqBW       = bytesNeeded / (ms * 1.0e-3) / 1.0e9;
    const double ai          = flops / bytesNeeded;             // = 0.25 FLOP/Byte
    const double roofCap     = peakBW * ai;
    const double peakGF      = peakFP32 * 1.0e3;

    printf("\n== Correctness ==\n");
    printf("  max abs error        : %.6e  (relative %.3e)\n", maxAbsErr, relErr);
    printf("  CPU reference time   : %.3f s\n", cpuSec);
    printf("  verification         : %s\n", (relErr < 1e-4) ? "PASS" : "FAIL");

    printf("\n== Performance (average of %d runs) ==\n", TIMED_ITERS);
    printf("  kernel time          : %.4f ms\n", ms);
    printf("  GFLOP/s              : %.2f\n", gflops);
    printf("  arithmetic intensity : %.3f FLOP/Byte  (4 Byte per FLOP)\n", ai);
    printf("  requested bytes      : %.3f GB  (2*N^3 loads x 4 Byte)\n",
           bytesNeeded / 1.0e9);
    printf("  requested bandwidth  : %.1f GB/s  (%.1f%% of peak)\n",
           reqBW, 100.0 * reqBW / peakBW);

    printf("\n== Roofline ==\n");
    printf("  machine balance      : %.1f FLOP/Byte vs naive AI %.2f -> memory bound\n",
           peakFP32 * 1e12 / (peakBW * 1e9), ai);
    printf("  roofline ceiling     : %.1f GFLOP/s  (%.1f GB/s x %.2f FLOP/Byte)\n",
           roofCap, peakBW, ai);
    printf("  ceiling as %% of peak: %.2f%%\n", 100.0 * roofCap / peakGF);
    printf("  measured / ceiling   : %.2f%%   <-- 实测性能占 Roofline 上限的比例\n",
           100.0 * gflops / roofCap);
    printf("  measured / peak      : %.2f%%\n", 100.0 * gflops / peakGF);

    // ---------- 9. 释放资源 ----------
    CUDA_CHECK(cudaEventDestroy(evStart));
    CUDA_CHECK(cudaEventDestroy(evStop));
    CUDA_CHECK(cudaFree(d_M));
    CUDA_CHECK(cudaFree(d_N));
    CUDA_CHECK(cudaFree(d_P));
    free(h_M); free(h_N); free(h_P); free(h_Ref);
    CUDA_CHECK(cudaDeviceReset());
    return 0;
}

// 文件: matmul_diag.cu
// 主题: ECE408 Lecture 4 -- 诊断实验：合并访问 vs 跨步访问（stride = Width）的带宽代价
// 编译: nvcc -O3 -arch=sm_80 matmul_diag.cu -o matmul_diag
// 运行: ./matmul_diag 4096
//       参数为探测矩阵边长 W，默认 4096（单矩阵 64 MB，远超 A100 的 40 MB L2，
//       只有把工作集赶出 cache 才能观察到真正的 DRAM 行为）

#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 PROBE_BLOCK   256
#define TILE_W         16
#define STREAM_ITERS   10

// ---------------- 探针 1：完全合并的顺序读（每 warp 128 Byte 连续） ----------------
__global__ void streamReadKernel(const float* __restrict__ in,
                                 float* __restrict__ out, int n)
{
    int g = blockIdx.x * blockDim.x + threadIdx.x;
    if (g < n) out[g] = in[g] * 2.0f;
}

// ---------------- 探针 2：跨步读（相邻线程地址相距 Width*4 Byte） ----------------
__global__ void stridedReadKernel(const float* __restrict__ in,
                                  float* __restrict__ out, int Width)
{
    const int g = blockIdx.x * blockDim.x + threadIdx.x;
    const int n = Width * Width;
    if (g < n) {
        const int Row = g / Width;    // 与 thread-per-output 中的输出行号同构
        const int Col = g % Width;    // 与 thread-per-output 中的输出列号同构
        out[g] = in[Col * Width + Row] * 2.0f;   // 相邻线程（Col+1）地址 + Width*4 Byte
    }
}

// ---------------- 探针 3：朴素矩阵乘法的访存流量（保留访问模式，去掉 FLOP） ----------------
__global__ void naiveLoadsOnlyKernel(const float* __restrict__ d_M,
                                     const float* __restrict__ d_N,
                                     float* __restrict__ out, int Width)
{
    const int Row = blockIdx.y * blockDim.y + threadIdx.y;
    const int Col = blockIdx.x * blockDim.x + threadIdx.x;
    if (Row < Width && Col < Width) {
        float acc = 0.0f;
        for (int k = 0; k < Width; ++k) {
            acc += d_M[Row * Width + k] + d_N[k * Width + Col];
        }
        out[Row * Width + Col] = acc;
    }
}

// ============================== 校验（抽样） ==============================
static void verifyStream(const float* h_out, const float* h_in, int n)
{
    int bad = 0;
    for (int g = 0; g < n; g += 9973) {
        if (h_out[g] != h_in[g] * 2.0f) ++bad;
    }
    printf("  verification         : %s (sampled)\n", (bad == 0) ? "PASS" : "FAIL");
}

static void verifyStrided(const float* h_out, const float* h_in, int W)
{
    const int n = W * W;
    int bad = 0;
    for (int g = 0; g < n; g += 9973) {
        const int Row = g / W;
        const int Col = g % W;
        if (h_out[g] != h_in[Col * W + Row] * 2.0f) ++bad;
    }
    printf("  verification         : %s (sampled)\n", (bad == 0) ? "PASS" : "FAIL");
}

static void verifyLoadsOnly(const float* h_out, const float* h_in, int W)
{
    const int n = W * W;
    int bad = 0;
    for (int g = 0; g < n; g += 262147) {
        const int Row = g / W;
        const int Col = g % W;
        double expect = 0.0;
        for (int k = 0; k < W; ++k) {
            expect += (double)h_in[Row * W + k] + (double)h_in[k * W + Col];
        }
        const double diff = (double)h_out[g] - expect;
        if (fabs(diff) > 1e-3 * fabs(expect) + 1.0) ++bad;
    }
    printf("  verification         : %s (sampled)\n", (bad == 0) ? "PASS" : "FAIL");
}

// ================================== main ===================================
int main(int argc, char** argv)
{
    const int    W     = (argc > 1) ? atoi(argv[1]) : 4096;
    const int    n     = W * W;
    const size_t bytes = (size_t)W * (size_t)W * sizeof(float);

    CUDA_CHECK(cudaSetDevice(0));
    cudaDeviceProp prop;
    CUDA_CHECK(cudaGetDeviceProperties(&prop, 0));
    const double peakBW = 2.0 * (double)prop.memoryClockRate * 1.0e3
                        * ((double)prop.memoryBusWidth / 8.0) / 1.0e9;

    printf("== Device ==\n");
    printf("  %s\n", prop.name);
    printf("  peak bandwidth       : %.1f GB/s\n", peakBW);
    printf("  L2 cache             : %.1f MB\n", prop.l2CacheSize / 1048576.0);
    printf("  probe width W        : %d   (one matrix = %.1f MB)\n",
           W, bytes / 1048576.0);
    printf("  stride between adjacent threads in probe 2: W*4 = %d Byte\n", W * 4);

    // ---------- 主机数据 ----------
    float* h_in  = (float*)malloc(bytes);
    float* h_out = (float*)malloc(bytes);
    if (!h_in || !h_out) { fprintf(stderr, "host malloc failed\n"); exit(EXIT_FAILURE); }
    for (int i = 0; i < n; ++i) h_in[i] = (float)(i % 17) * 0.125f;

    float *d_in = NULL, *d_out = NULL;
    CUDA_CHECK(cudaMalloc((void**)&d_in, bytes));
    CUDA_CHECK(cudaMalloc((void**)&d_out, bytes));
    CUDA_CHECK(cudaMemcpy(d_in, h_in, bytes, cudaMemcpyHostToDevice));

    cudaEvent_t evA, evB;
    CUDA_CHECK(cudaEventCreate(&evA));
    CUDA_CHECK(cudaEventCreate(&evB));

    const int gridStream = (n + PROBE_BLOCK - 1) / PROBE_BLOCK;
    const double usefulBytes = (double)bytes;      // 每个探针"有用"的读字节数

    // ================= 探针 1：顺序（合并）读 =================
    for (int i = 0; i < 2; ++i) {
        streamReadKernel<<<gridStream, PROBE_BLOCK>>>(d_in, d_out, n);
    }
    CUDA_CHECK(cudaDeviceSynchronize());
    CUDA_CHECK(cudaEventRecord(evA));
    for (int i = 0; i < STREAM_ITERS; ++i) {
        streamReadKernel<<<gridStream, PROBE_BLOCK>>>(d_in, d_out, n);
    }
    CUDA_CHECK(cudaEventRecord(evB));
    CUDA_CHECK(cudaEventSynchronize(evB));
    float ms1 = 0.0f;
    CUDA_CHECK(cudaEventElapsedTime(&ms1, evA, evB));
    ms1 /= STREAM_ITERS;
    CUDA_CHECK(cudaMemcpy(h_out, d_out, bytes, cudaMemcpyDeviceToHost));
    const double bw1 = usefulBytes / (ms1 * 1.0e-3) / 1.0e9;

    printf("\n== Probe 1: coalesced sequential read ==\n");
    printf("  read bytes per iter  : %.1f MB\n", usefulBytes / 1048576.0);
    printf("  time                 : %.4f ms\n", ms1);
    printf("  achieved bandwidth   : %.1f GB/s  (%.1f%% of peak)\n",
           bw1, 100.0 * bw1 / peakBW);
    printf("  model                : 32 threads x 4 B = 128 B = 4 sectors, 100%% useful\n");
    verifyStream(h_out, h_in, n);

    // ================= 探针 2：跨步读 =================
    for (int i = 0; i < 2; ++i) {
        stridedReadKernel<<<gridStream, PROBE_BLOCK>>>(d_in, d_out, W);
    }
    CUDA_CHECK(cudaDeviceSynchronize());
    CUDA_CHECK(cudaEventRecord(evA));
    for (int i = 0; i < STREAM_ITERS; ++i) {
        stridedReadKernel<<<gridStream, PROBE_BLOCK>>>(d_in, d_out, W);
    }
    CUDA_CHECK(cudaEventRecord(evB));
    CUDA_CHECK(cudaEventSynchronize(evB));
    float ms2 = 0.0f;
    CUDA_CHECK(cudaEventElapsedTime(&ms2, evA, evB));
    ms2 /= STREAM_ITERS;
    CUDA_CHECK(cudaMemcpy(h_out, d_out, bytes, cudaMemcpyDeviceToHost));
    const double bw2 = usefulBytes / (ms2 * 1.0e-3) / 1.0e9;

    printf("\n== Probe 2: strided read (stride = W*4 Byte) ==\n");
    printf("  read bytes per iter  : %.1f MB\n", usefulBytes / 1048576.0);
    printf("  time                 : %.4f ms\n", ms2);
    printf("  achieved bandwidth   : %.1f GB/s  (%.1f%% of peak)\n",
           bw2, 100.0 * bw2 / peakBW);
    printf("  model                : 32 addresses spaced %d B apart\n", W * 4);
    printf("                         -> 32 sectors x 32 B = 1024 B fetched, 128 B useful\n");
    printf("                         -> sector efficiency 12.5%%, i.e. 8x waste\n");
    printf("  slowdown vs probe 1  : %.2fx\n", ms2 / ms1);
    verifyStrided(h_out, h_in, W);

    // ================= 探针 3：朴素矩阵乘法的访存流量 =================
    dim3 dimBlock(TILE_W, TILE_W, 1);
    dim3 dimGrid((W + TILE_W - 1) / TILE_W, (W + TILE_W - 1) / TILE_W, 1);
    CUDA_CHECK(cudaEventRecord(evA));
    naiveLoadsOnlyKernel<<<dimGrid, dimBlock>>>(d_in, d_in, d_out, W);
    CUDA_CHECK(cudaEventRecord(evB));
    CUDA_CHECK(cudaEventSynchronize(evB));
    float ms3 = 0.0f;
    CUDA_CHECK(cudaEventElapsedTime(&ms3, evA, evB));
    CUDA_CHECK(cudaMemcpy(h_out, d_out, bytes, cudaMemcpyDeviceToHost));

    const double requested3  = 2.0 * (double)W * W * W * 4.0;      // 2*W^3 次 load
    const double warps3      = (double)(dimGrid.x * dimGrid.y) * 8.0;
    const double modelFetch3 = warps3 * (double)W * 128.0;         // 每 warp 每 k 128 B
    const double bw3         = requested3 / (ms3 * 1.0e-3) / 1.0e9;

    printf("\n== Probe 3: naive matmul access pattern without FMA ==\n");
    printf("  grid                 : (%u, %u, 1), block (%u, %u, 1)\n",
           dimGrid.x, dimGrid.y, dimBlock.x, dimBlock.y);
    printf("  requested bytes      : %.2f GB  (2*W^3 loads x 4 Byte)\n",
           requested3 / 1.0e9);
    printf("  modeled fetched bytes: %.2f GB  (warps x W x 128 B per warp per k)\n",
           modelFetch3 / 1.0e9);
    printf("  time                 : %.4f ms\n", ms3);
    printf("  requested bandwidth  : %.1f GB/s  (%.1f%% of peak)\n",
           bw3, 100.0 * bw3 / peakBW);
    printf("  DRAM-bound lower bnd : %.4f ms  (modeled fetch / peak BW)\n",
           modelFetch3 / (peakBW * 1e9) * 1.0e3);
    verifyLoadsOnly(h_out, h_in, W);

    // ---------- 释放 ----------
    CUDA_CHECK(cudaEventDestroy(evA));
    CUDA_CHECK(cudaEventDestroy(evB));
    CUDA_CHECK(cudaFree(d_in));
    CUDA_CHECK(cudaFree(d_out));
    free(h_in);
    free(h_out);
    CUDA_CHECK(cudaDeviceReset());
    return 0;
}

// 文件: matmul_transposed.cu
// 主题: ECE408 Lecture 4 -- 存储布局对合并访存的影响（A 转置 vs B 转置）
// 编译: nvcc -O3 -arch=sm_80 matmul_transposed.cu -o matmul_transposed
// 运行: ./matmul_transposed 2048
//       参数为矩阵边长 Width（默认 2048，单矩阵 16 MB，三个矩阵共 48 MB > A100 L2 40 MB）

#include <cstdio>
#include <cstdlib>
#include <cmath>
#include <ctime>
#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 BLOCK_WIDTH  16
#define WARMUP_ITERS  3
#define TIMED_ITERS  10

// ============================== device kernels =============================
// (a) 转置 kernel：out[c*Width + r] = in[r*Width + c]
//     读：相邻 tx -> 相邻地址（合并）；写：相邻 tx 的地址相差 Width*4 Byte（跨步）
__global__ void transposeKernel(const float* __restrict__ in,
                                float* __restrict__ out, int Width)
{
    const int Row = blockIdx.y * blockDim.y + threadIdx.y;
    const int Col = blockIdx.x * blockDim.x + threadIdx.x;
    if (Row < Width && Col < Width) {
        out[Col * Width + Row] = in[Row * Width + Col];
    }
}

// (b) 朴素版：A、B 都是行主序。A 侧跨步（2 个 sector），B 侧合并（2 个 sector）
__global__ void matmulNaiveKernel(const float* __restrict__ d_M,
                                  const float* __restrict__ d_N,
                                  float* __restrict__ d_P, int Width)
{
    const int Row = blockIdx.y * blockDim.y + threadIdx.y;
    const int Col = blockIdx.x * blockDim.x + threadIdx.x;
    if (Row < Width && Col < Width) {
        float Pvalue = 0.0f;
        for (int k = 0; k < Width; ++k) {
            Pvalue += d_M[Row * Width + k] * d_N[k * Width + Col];
        }
        d_P[Row * Width + Col] = Pvalue;
    }
}

// (c) A 转置版：M_T[k][Row] = M[Row][k]。A 侧变成"2 个相邻地址"（1 个 sector）
__global__ void matmulAMTKernel(const float* __restrict__ d_MT,
                                const float* __restrict__ d_N,
                                float* __restrict__ d_P, int Width)
{
    const int Row = blockIdx.y * blockDim.y + threadIdx.y;
    const int Col = blockIdx.x * blockDim.x + threadIdx.x;
    if (Row < Width && Col < Width) {
        float Pvalue = 0.0f;
        for (int k = 0; k < Width; ++k) {
            Pvalue += d_MT[k * Width + Row] * d_N[k * Width + Col];
        }
        d_P[Row * Width + Col] = Pvalue;
    }
}

// (d) B 转置版：N_T[Col][k] = N[k][Col]，即把 B 按列主序存储后按行读取
//     B 侧变成"16 个相距 Width*4 Byte 的地址"（16 个独立 sector）
__global__ void matmulBTKernel(const float* __restrict__ d_M,
                               const float* __restrict__ d_NT,
                               float* __restrict__ d_P, int Width)
{
    const int Row = blockIdx.y * blockDim.y + threadIdx.y;
    const int Col = blockIdx.x * blockDim.x + threadIdx.x;
    if (Row < Width && Col < Width) {
        float Pvalue = 0.0f;
        for (int k = 0; k < Width; ++k) {
            Pvalue += d_M[Row * Width + k] * d_NT[Col * Width + k];
        }
        d_P[Row * Width + Col] = Pvalue;
    }
}

// ============================== host utilities =============================
static void initMatrix(float* a, int n, unsigned int seed)
{
    srand(seed);
    for (int i = 0; i < n * n; ++i) a[i] = (float)(rand() % 7) * 0.25f;
}

static void cpuMatmul(const float* A, const float* B, float* C, int n)
{
    for (int i = 0; i < n * n; ++i) C[i] = 0.0f;
    for (int i = 0; i < n; ++i) {
        for (int k = 0; k < n; ++k) {
            const float  a    = A[i * n + k];
            const float* Brow = B + k * n;
            float*       Crow = C + i * n;
            for (int j = 0; j < n; ++j) Crow[j] += a * Brow[j];
        }
    }
}

static void cpuTranspose(const float* in, float* out, int n)
{
    for (int r = 0; r < n; ++r) {
        for (int c = 0; c < n; ++c) {
            out[c * n + r] = in[r * n + c];
        }
    }
}

static int fp32CoresPerSM(int major, int minor)
{
    if (major == 7) return 64;
    if (major == 8) return (minor == 0) ? 64 : 128;
    if (major >= 9) return 128;
    return 128;
}

// 计时辅助：把一次 kernel 启动包成可调用对象，返回平均毫秒数
template <typename LaunchFn>
static double timeKernel(LaunchFn launch, int warmup, int iters)
{
    for (int i = 0; i < warmup; ++i) launch();
    CUDA_CHECK(cudaGetLastError());
    CUDA_CHECK(cudaDeviceSynchronize());
    cudaEvent_t a, b;
    CUDA_CHECK(cudaEventCreate(&a));
    CUDA_CHECK(cudaEventCreate(&b));
    CUDA_CHECK(cudaEventRecord(a));
    for (int i = 0; i < iters; ++i) launch();
    CUDA_CHECK(cudaEventRecord(b));
    CUDA_CHECK(cudaEventSynchronize(b));
    float ms = 0.0f;
    CUDA_CHECK(cudaEventElapsedTime(&ms, a, b));
    CUDA_CHECK(cudaEventDestroy(a));
    CUDA_CHECK(cudaEventDestroy(b));
    return (double)ms / iters;
}

static void verifyMatrix(const float* got, const float* want, int n, const char* tag)
{
    double maxAbs = 0.0, maxRef = 0.0;
    for (int i = 0; i < n * n; ++i) {
        const double e = fabs((double)got[i] - (double)want[i]);
        if (e > maxAbs) maxAbs = e;
        const double r = fabs((double)want[i]);
        if (r > maxRef) maxRef = r;
    }
    const double rel = (maxRef > 0.0) ? maxAbs / maxRef : maxAbs;
    printf("  %-26s: max abs err %.3e, rel %.3e -> %s\n",
           tag, maxAbs, rel, (rel < 1e-4) ? "PASS" : "FAIL");
}

// ================================== main ===================================
int main(int argc, char** argv)
{
    const int    N     = (argc > 1) ? atoi(argv[1]) : 2048;
    const size_t bytes = (size_t)N * (size_t)N * sizeof(float);

    CUDA_CHECK(cudaSetDevice(0));
    cudaDeviceProp prop;
    CUDA_CHECK(cudaGetDeviceProperties(&prop, 0));
    const double peakBW   = 2.0 * (double)prop.memoryClockRate * 1.0e3
                          * ((double)prop.memoryBusWidth / 8.0) / 1.0e9;
    const double peakFP32 = (double)prop.multiProcessorCount
                          * (double)fp32CoresPerSM(prop.major, prop.minor)
                          * 2.0 * (double)prop.clockRate * 1.0e3 / 1.0e12;
    const double flops    = 2.0 * (double)N * (double)N * (double)N;

    printf("== Device ==\n");
    printf("  %s, peak BW %.1f GB/s, peak FP32 %.2f TFLOPS, L2 %.1f MB\n",
           prop.name, peakBW, peakFP32, prop.l2CacheSize / 1048576.0);
    printf("  N = %d, three matrices = %.1f MB\n", N, 3.0 * bytes / 1048576.0);

    // ---------- 主机端数据与参考结果 ----------
    float* h_M   = (float*)malloc(bytes);
    float* h_N   = (float*)malloc(bytes);
    float* h_MT  = (float*)malloc(bytes);
    float* h_NT  = (float*)malloc(bytes);
    float* h_P   = (float*)malloc(bytes);
    float* h_Ref = (float*)malloc(bytes);
    if (!h_M || !h_N || !h_MT || !h_NT || !h_P || !h_Ref) {
        fprintf(stderr, "host malloc failed\n");
        exit(EXIT_FAILURE);
    }
    initMatrix(h_M, N, 3u);
    initMatrix(h_N, N, 4u);
    cpuTranspose(h_M, h_MT, N);
    cpuTranspose(h_N, h_NT, N);

    const clock_t t0 = clock();
    cpuMatmul(h_M, h_N, h_Ref, N);
    const double cpuSec = (double)(clock() - t0) / CLOCKS_PER_SEC;

    // ---------- 设备端内存 ----------
    float *d_M = NULL, *d_N = NULL, *d_MT = NULL, *d_NT = NULL, *d_P = NULL;
    CUDA_CHECK(cudaMalloc((void**)&d_M,  bytes));
    CUDA_CHECK(cudaMalloc((void**)&d_N,  bytes));
    CUDA_CHECK(cudaMalloc((void**)&d_MT, bytes));
    CUDA_CHECK(cudaMalloc((void**)&d_NT, bytes));
    CUDA_CHECK(cudaMalloc((void**)&d_P,  bytes));
    CUDA_CHECK(cudaMemcpy(d_M, h_M, bytes, cudaMemcpyHostToDevice));
    CUDA_CHECK(cudaMemcpy(d_N, h_N, bytes, cudaMemcpyHostToDevice));

    dim3 dimBlock(BLOCK_WIDTH, BLOCK_WIDTH, 1);
    dim3 dimGrid((N + BLOCK_WIDTH - 1) / BLOCK_WIDTH,
                 (N + BLOCK_WIDTH - 1) / BLOCK_WIDTH, 1);

    // ---------- 转置（GPU） ----------
    const double msT_M = timeKernel([&] {
        transposeKernel<<<dimGrid, dimBlock>>>(d_M, d_MT, N);
    }, WARMUP_ITERS, TIMED_ITERS);
    const double msT_N = timeKernel([&] {
        transposeKernel<<<dimGrid, dimBlock>>>(d_N, d_NT, N);
    }, WARMUP_ITERS, TIMED_ITERS);
    const double transposeGBs = 2.0 * (double)bytes / (msT_M * 1.0e-3) / 1.0e9;

    CUDA_CHECK(cudaMemcpy(h_P, d_MT, bytes, cudaMemcpyDeviceToHost));
    printf("\n== Transpose kernel ==\n");
    printf("  time per transpose   : %.4f ms (M_T), %.4f ms (N_T)\n", msT_M, msT_N);
    printf("  effective bandwidth  : %.1f GB/s (read+write counted)\n", transposeGBs);
    printf("  model                : read coalesced (2 sectors), write strided (16 sectors)\n");
    verifyMatrix(h_P, h_MT, N, "M_T = transpose(M)");
    CUDA_CHECK(cudaMemcpy(h_P, d_NT, bytes, cudaMemcpyDeviceToHost));
    verifyMatrix(h_P, h_NT, N, "N_T = transpose(N)");

    // ---------- 三个矩阵乘法 kernel ----------
    const double msNaive = timeKernel([&] {
        matmulNaiveKernel<<<dimGrid, dimBlock>>>(d_M, d_N, d_P, N);
    }, WARMUP_ITERS, TIMED_ITERS);
    CUDA_CHECK(cudaMemcpy(h_P, d_P, bytes, cudaMemcpyDeviceToHost));
    printf("\n== matmulNaiveKernel (A row-major, B row-major) ==\n");
    printf("  time                 : %.4f ms   ->  %.2f GFLOP/s\n",
           msNaive, flops / (msNaive * 1.0e-3) / 1.0e9);
    printf("  non-coalesced side   : A  (2 sectors = 64 B fetched for 8 B useful)\n");
    verifyMatrix(h_P, h_Ref, N, "P = M*N (naive)");

    const double msAMT = timeKernel([&] {
        matmulAMTKernel<<<dimGrid, dimBlock>>>(d_MT, d_N, d_P, N);
    }, WARMUP_ITERS, TIMED_ITERS);
    CUDA_CHECK(cudaMemcpy(h_P, d_P, bytes, cudaMemcpyDeviceToHost));
    printf("\n== matmulAMTKernel (A transposed) ==\n");
    printf("  time                 : %.4f ms   ->  %.2f GFLOP/s\n",
           msAMT, flops / (msAMT * 1.0e-3) / 1.0e9);
    printf("  A side               : 2 adjacent addresses -> 1 sector = 32 B (was 64 B)\n");
    verifyMatrix(h_P, h_Ref, N, "P = M*N (A transposed)");

    const double msBT = timeKernel([&] {
        matmulBTKernel<<<dimGrid, dimBlock>>>(d_M, d_NT, d_P, N);
    }, WARMUP_ITERS, TIMED_ITERS);
    CUDA_CHECK(cudaMemcpy(h_P, d_P, bytes, cudaMemcpyDeviceToHost));
    printf("\n== matmulBTKernel (B stored column-major) ==\n");
    printf("  time                 : %.4f ms   ->  %.2f GFLOP/s\n",
           msBT, flops / (msBT * 1.0e-3) / 1.0e9);
    printf("  B side               : 16 addresses spaced Width*4 B -> 16 sectors = 512 B\n");
    verifyMatrix(h_P, h_Ref, N, "P = M*N (B column-major)");

    // ---------- 汇总 ----------
    printf("\n== Summary (N = %d) ==\n", N);
    printf("  kernel                        ms      GFLOP/s   warp-bytes per k   sector eff\n");
    printf("  naive   (A row, B row)     %7.4f   %8.2f        128 B               56%%\n",
           msNaive, flops / (msNaive * 1.0e-3) / 1.0e9);
    printf("  A-trans (A^T, B row)       %7.4f   %8.2f         96 B               75%%\n",
           msAMT, flops / (msAMT * 1.0e-3) / 1.0e9);
    printf("  B-trans (A row, B col)     %7.4f   %8.2f        576 B               12.5%%\n",
           msBT, flops / (msBT * 1.0e-3) / 1.0e9);
    printf("  ratio B-trans / A-trans    %.2fx\n", msBT / msAMT);
    printf("  ratio B-trans / naive      %.2fx\n", msBT / msNaive);
    printf("  ratio naive   / A-trans    %.2fx\n", msNaive / msAMT);

    printf("\n== Roofline view (naive) ==\n");
    printf("  arithmetic intensity : 0.25 FLOP/Byte (4 Byte per FLOP)\n");
    printf("  roofline ceiling     : %.1f GFLOP/s (%.1f GB/s x 0.25)\n",
           peakBW * 0.25, peakBW);
    printf("  ceiling / peak FP32  : %.2f%%\n", 100.0 * peakBW * 0.25 / (peakFP32 * 1.0e3));
    printf("  CPU reference time   : %.3f s\n", cpuSec);

    // ---------- 释放 ----------
    CUDA_CHECK(cudaFree(d_M));
    CUDA_CHECK(cudaFree(d_N));
    CUDA_CHECK(cudaFree(d_MT));
    CUDA_CHECK(cudaFree(d_NT));
    CUDA_CHECK(cudaFree(d_P));
    free(h_M); free(h_N); free(h_MT); free(h_NT); free(h_P); free(h_Ref);
    CUDA_CHECK(cudaDeviceReset());
    return 0;
}