📚 课程目录第2章矩阵乘法基础实现

2. 矩阵乘法基础实现

任务目标

编写一个基础的GPU矩阵乘法内核（C = A × B），作为后续优化的性能基线。内核不包含任何优化技巧（无共享内存、无向量化、无循环展开），仅实现最直接的并行策略：每个线程计算输出矩阵中的一个元素。同时提供主机端代码完成数据初始化、内核启动、结果校验与耗时测量。

环境与依赖

CUDA Toolkit 11.0+（支持 cudaEvent 计时）
GPU 计算能力 6.0+（兼容任意现代 NVIDIA GPU）
编译命令：nvcc -O2 -arch=sm_70 matmul_baseline.cu -o matmul_baseline

关键完整代码

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

#define CHECK_CUDA(call) do { \
    cudaError_t err = (call); \
    if (err != cudaSuccess) { \
        printf("CUDA error at %s:%d: %s\n", __FILE__, __LINE__, \
               cudaGetErrorString(err)); \
        exit(1); \
    } \
} while(0)

// 基线内核：每个线程计算 C 的一个元素
__global__ void matmul_naive(const float* A, const float* B, float* C,
                             int M, int N, int K) {
    int row = blockIdx.y * blockDim.y + threadIdx.y;
    int col = blockIdx.x * blockDim.x + threadIdx.x;
    if (row < M && col < N) {
        float sum = 0.0f;
        for (int k = 0; k < K; ++k) {
            sum += A[row * K + k] * B[k * N + col];
        }
        C[row * N + col] = sum;
    }
}

int main() {
    const int M = 1024, N = 1024, K = 1024;
    size_t bytes_A = M * K * sizeof(float);
    size_t bytes_B = K * N * sizeof(float);
    size_t bytes_C = M * N * sizeof(float);

    // 主机端分配与初始化
    float *h_A = (float*)malloc(bytes_A);
    float *h_B = (float*)malloc(bytes_B);
    float *h_C = (float*)malloc(bytes_C);
    float *h_C_ref = (float*)malloc(bytes_C);
    for (int i = 0; i < M * K; ++i) h_A[i] = 1.0f;
    for (int i = 0; i < K * N; ++i) h_B[i] = 0.001f;
    for (int i = 0; i < M * N; ++i) h_C_ref[i] = K * 0.001f; // 预期结果

    // 设备端分配
    float *d_A, *d_B, *d_C;
    CHECK_CUDA(cudaMalloc(&d_A, bytes_A));
    CHECK_CUDA(cudaMalloc(&d_B, bytes_B));
    CHECK_CUDA(cudaMalloc(&d_C, bytes_C));
    CHECK_CUDA(cudaMemcpy(d_A, h_A, bytes_A, cudaMemcpyHostToDevice));
    CHECK_CUDA(cudaMemcpy(d_B, h_B, bytes_B, cudaMemcpyHostToDevice));

    // 内核启动配置：16x16 线程块
    dim3 block(16, 16);
    dim3 grid((N + block.x - 1) / block.x, (M + block.y - 1) / block.y);

    // 预热 + 计时
    matmul_naive<<<grid, block>>>(d_A, d_B, d_C, M, N, K);
    CHECK_CUDA(cudaDeviceSynchronize());

    cudaEvent_t start, stop;
    CHECK_CUDA(cudaEventCreate(&start));
    CHECK_CUDA(cudaEventCreate(&stop));
    CHECK_CUDA(cudaEventRecord(start));
    matmul_naive<<<grid, block>>>(d_A, d_B, d_C, M, N, K);
    CHECK_CUDA(cudaEventRecord(stop));
    CHECK_CUDA(cudaEventSynchronize(stop));

    float ms = 0.0f;
    CHECK_CUDA(cudaEventElapsedTime(&ms, start, stop));
    printf("Kernel time: %.3f ms\n", ms);

    // 结果校验
    CHECK_CUDA(cudaMemcpy(h_C, d_C, bytes_C, cudaMemcpyDeviceToHost));
    float max_err = 0.0f;
    for (int i = 0; i < M * N; ++i) {
        float err = fabsf(h_C[i] - h_C_ref[i]);
        if (err > max_err) max_err = err;
    }
    printf("Max error: %e\n", max_err);
    printf(max_err < 1e-4 ? "PASS\n" : "FAIL\n");

    // 清理
    free(h_A); free(h_B); free(h_C); free(h_C_ref);
    cudaFree(d_A); cudaFree(d_B); cudaFree(d_C);
    cudaEventDestroy(start); cudaEventDestroy(stop);
    return 0;
}

运行与验证

保存代码为 matmul_baseline.cu，按上述命令编译。
运行 ./matmul_baseline，观察输出。
通过条件：程序输出 PASS，且打印的 Kernel time 为正值（具体数值取决于 GPU 型号，通常为几十毫秒量级）。

验证逻辑说明

输入数据采用全 1 矩阵 A 与全 0.001 矩阵 B，因此 C 的每个元素理论值为 K × 0.001 = 1.024。校验时比较 GPU 计算结果与理论值的最大误差，误差小于 1e-4 即判定正确。该设计使校验不依赖额外 CPU 矩阵乘法实现，同时能精确捕获计算错误。

上图展示了基线内核的并行映射关系：每个线程负责输出矩阵 C 中一个固定位置的元素，通过遍历 K 维度完成点积累加。该实现无数据复用，每个 A 元素被读取 N 次、每个 B 元素被读取 M 次，访存效率低，这正是后续优化（如共享内存分块）的切入点。