2. CUDA编程模型与矩阵运算基础

CUDA 编程模型的核心思想，是把一个大规模数据并行任务拆解成大量可以同时执行的小任务，并把这些任务映射到 GPU 的硬件执行单元上。理解这一模型，需要同时把握三条主线：线程如何组织、内存如何分层、核函数如何编写。这三者共同决定了矩阵运算在 GPU 上的执行效率。

线程层次：从网格到线程

CUDA 用三级结构组织线程：网格（Grid）、线程块（Block）、线程（Thread）。一个核函数启动时对应一个网格，网格由若干线程块组成，每个线程块又包含若干线程。线程块内的线程可以通过共享内存和同步指令协作，而不同线程块之间默认无法直接通信。

每个线程拥有唯一的三维索引：blockIdx 表示所在线程块在网格中的位置，threadIdx 表示线程在线程块内的位置，blockDim 表示每个线程块的维度。对于一维情况，全局线程编号常用 blockIdx.x * blockDim.x + threadIdx.x 计算。这个编号就是线程处理数据的“身份证”，矩阵元素的访问几乎都依赖它。

线程块大小需要权衡：太小则并行度不足、调度开销占比高；太大则受寄存器与共享内存限制，可能降低占用率。常见做法是取 128、256 或 512 个线程，并尽量让线程块大小是 32（一个 warp）的整数倍，以避免最后一个 warp 中大量线程闲置。

内存层次：决定性能的关键

GPU 的内存层次直接对应线程层次。每个线程拥有私有的寄存器，速度最快但数量有限；每个线程块拥有一块共享内存，块内线程可读写，常用于缓存被反复使用的数据；所有线程共享全局内存，容量大但延迟高，是主机与设备之间数据传输的主要载体。此外还有只读的常量内存和纹理内存，适合存放核函数参数或查表数据。

矩阵运算的性能瓶颈往往不在计算，而在全局内存访问。以矩阵加法为例，每个输出元素只读取两个输入元素、写入一个元素，计算量极小，属于典型的内存带宽受限问题。此时应保证相邻线程访问相邻的全局内存地址，即合并访问，让一个 warp 的 32 次访问合并成少量内存事务。若线程按列优先方式跨步访问行优先存储的矩阵，带宽利用率会显著下降。

矩阵乘法则是另一类问题。朴素实现中，每个输出元素需要读取一整行和一整列，全局内存访问次数是 O(N³)。通过把矩阵分块载入共享内存，让块内线程复用这些数据，可以把全局内存访问降到 O(N³/B) 量级，B 为分块边长。这正是共享内存在矩阵乘法中的核心价值。

核函数编写：从加法到乘法

核函数用 __global__ 修饰，由主机端通过 <<<grid, block>>> 语法启动，在设备端由多个线程并行执行。下面是一个矩阵加法的核函数示例，矩阵按行优先存储，尺寸为 M×N：

__global__ void matAdd(const float* A, const float* B, float* C, int M, int N) {
    int col = blockIdx.x * blockDim.x + threadIdx.x;
    int row = blockIdx.y * blockDim.y + threadIdx.y;
    if (row < M && col < N) {
        int idx = row * N + col;
        C[idx] = A[idx] + B[idx];
    }
}

主机端启动时，通常把线程块设为 dim3 block(16, 16)，网格设为 dim3 grid((N + 15) / 16, (M + 15) / 16)。边界判断 if (row < M && col < N) 不可省略，因为网格尺寸向上取整后，最后一列或最后一行会多出无效线程。以 1000×1000 矩阵为例，block 为 16×16 时，grid 为 63×63，实际覆盖 1008×1008 个线程，多出的 8 行 8 列必须被条件拦截。

矩阵乘法的朴素核函数中，每个线程计算 C 的一个元素，沿 K 维循环累加：

__global__ void matMulNaive(const float* A, const float* B, float* C,
                            int M, int N, int K) {
    int col = blockIdx.x * blockDim.x + threadIdx.x;
    int row = blockIdx.y * blockDim.y + threadIdx.y;
    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;
    }
}

这段代码正确但效率不高：内层循环中 A 的元素被同一行的多个线程重复读取，B 的元素被同一列的多个线程重复读取，全局内存流量远大于必要值。改进方向是把 A、B 的分块载入共享内存，每个线程只从全局内存读取一次分块数据，然后在共享内存上完成 B×B 次乘加。分块大小常取 16×16 或 32×32，需保证 blockDim 与分块尺寸一致，并在载入后调用 __syncthreads() 确保数据就绪。

主机与设备的数据流转

一个完整的矩阵运算流程包括：主机分配并初始化数据、在设备上分配显存（cudaMalloc）、把数据从主机拷贝到设备（cudaMemcpy 的 HostToDevice 方向）、启动核函数、把结果拷回主机（DeviceToHost 方向）、释放显存。核函数启动是异步的，若紧接着要读取结果，需用 cudaDeviceSynchronize() 或依赖拷贝操作的隐式同步。实际调试时，应在每次 CUDA API 调用后检查返回值，并用 cudaGetLastError() 捕获核函数启动配置错误，例如线程块超过 1024 个线程这类问题。

验证结果时，对矩阵加法可抽取若干元素与 CPU 计算结果比对；对矩阵乘法，由于浮点累加顺序不同，GPU 结果与 CPU 结果可能存在微小差异，通常用相对误差阈值（如 1e-5）判断，而不是要求逐位相等。以 512×512 单精度矩阵乘法为例，若 CPU 参考实现得到某元素为 128.375，GPU 结果在 128.3749 至 128.3751 之间即可认为通过。掌握线程映射、内存层次与核函数结构之后，后续的批处理与性能优化才有可靠的落脚点。