CUDA内存优化入门：从Naive到Tiled矩阵乘法—Week2学习总结

CUDA内存优化入门：从Naive到Tiled矩阵乘法—Week2学习总结

实战演示如何使用Shared Memory实现10倍+性能提升

一、为什么要关注矩阵乘法的优化？

如果你开始学习GPU编程，矩阵乘法几乎是一个绕不开的话题。这不仅仅因为它是一个经典的并行计算问题，更重要的是，矩阵乘法是深度学习的核心运算。无论是神经网络的前向传播、反向传播，还是注意力机制的计算，背后都离不开大量的矩阵乘法运算。

在现代深度学习框架中，矩阵乘法往往占据了70%以上的计算时间。这意味着，如果你能将矩阵乘法的性能提升一倍，整个神经网络的训练速度就能提升接近一倍。更令人兴奋的是，通过合理的优化，矩阵乘法的性能提升空间可以达到100倍甚至更多！

本文你将学到什么

通过这篇文章，你将：

• 理解GPU内存层级结构：了解为什么有些内存快，有些内存慢
• 掌握Shared Memory的使用：学会GPU编程中最重要的优化技术之一
• 实现三个版本的矩阵乘法：从CPU串行到GPU并行，再到GPU优化
• 看到真实的性能数据：亲眼见证10倍以上的性能提升

前置知识要求

在开始之前，你需要具备：

• Week 1的CUDA基础知识（Grid、Block、Thread的概念）
• 基本的C++编程能力
• 矩阵乘法的数学原理

如果你还不太熟悉CUDA的基本概念，建议先回顾一下Week 1的内容。

二、GPU内存层级：速度与容量的权衡

在开始优化之前，我们必须深入理解GPU的内存结构。这就像你在城市规划中需要了解不同的交通工具一样——飞机快但容量小，火车慢但能运更多货物。

三层内存结构

GPU的内存可以简单分为三个层次：

┌─────────────────────────────────────┐
│   Registers (寄存器)                │  最快，最小，每个线程私有
├─────────────────────────────────────┤
│   Shared Memory (共享内存)          │  快，小，block内共享 ← 本文重点！
├─────────────────────────────────────┤
│   Global Memory (全局内存)          │  慢，大，所有threads可访问
└─────────────────────────────────────┘

这三层内存就像你工作时的不同存储位置：

• Registers（寄存器）：就像你的大脑记忆，访问最快，但容量极小
• Shared Memory（共享内存）：就像你办公桌上的文件，团队成员都能快速访问
• Global Memory（全局内存）：就像公司的档案室，容量很大但需要走过去拿

性能对比数据

让我们看看这三种内存的具体性能差异：

从表格中可以看出，Shared Memory的访问速度比Global Memory快200倍左右！这就是我们优化的关键所在。

关键洞察：数据复用

为什么Shared Memory如此重要？核心在于数据复用。想象一下这个场景：

你需要从档案室（Global Memory）取一份文件，如果每次用都要跑过去拿，效率很低。更聪明的做法是：把这份文件复印一份放在办公桌上（Shared Memory），然后多次使用。这样，你只需要跑一次档案室，后续的访问都可以在办公桌上快速完成。

这就是我们要实现的优化思路：减少对Global Memory的访问，增加对Shared Memory的使用。

三、矩阵乘法的计算原理

在开始编写代码之前，让我们回顾一下矩阵乘法的基本原理。

数学公式

对于两个矩阵 A (M×K) 和 B (K×N)，它们的乘积 C (M×N) 的每个元素计算公式为：

C[i][j] = Σ(k=0 to K-1) A[i][k] * B[k][j]

用简单的话说：C矩阵第i行第j列的元素，等于A矩阵第i行与B矩阵第j列对应元素乘积的和。

计算量分析

对于规模为 N×N 的方阵乘法：

• 总乘法次数：N³ 次（每个元素需要N次乘法，共N²个元素）
• 总加法次数：N³ 次（每个元素需要N-1次加法，约等于N次）
• 计算复杂度：O(N³)

这意味着，当矩阵规模翻倍时，计算量会增加8倍！

内存访问分析：优化的关键

这是最重要的部分。让我们仔细分析内存访问模式：

计算C的一个元素需要：

• 读取A的一行（N个元素）
• 读取B的一列（N个元素）
• 总共：2N次内存读取

计算整个C矩阵需要：

• C有N²个元素
• 每个元素需要2N次读取
• 总共：2N³ 次内存读取

问题来了：A矩阵的每一行会被读取N次（因为C有N列），B矩阵的每一列会被读取N次（因为C有N行）。也就是说，同样的数据被重复读取了N次！

这就是优化的机会所在。如果我们能让数据在Shared Memory中被复用，就能大幅减少对Global Memory的访问次数。

三种实现方式的演进

• CPU版本：3层嵌套循环，简单直接，但完全串行
• GPU Naive版本：并行计算每个元素，但每次都从Global Memory读取
• GPU Tiled版本：使用Shared Memory缓存数据块，实现数据复用

接下来，我们将依次实现这三个版本，并对比它们的性能。

四、版本1：CPU基准实现

首先，我们实现一个最简单的CPU版本，作为性能对比的基准。

代码实现

void matmul_cpu(const float* A, const float* B, float* C, 
                int M, int N, int K) {
    // 外层循环：遍历C的每一行
    for (int i = 0; i < M; i++) {
        // 中层循环：遍历C的每一列
        for (int j = 0; j < N; j++) {
            float sum = 0.0f;
            // 内层循环：计算点积
            for (int k = 0; k < K; k++) {
                sum += A[i * K + k] * B[k * N + j];
            }
            C[i * N + j] = sum;
        }
    }
}

代码解析

这是最直观的矩阵乘法实现：

1. 外层两个循环确定C矩阵中要计算的元素位置(i, j)
2. 内层循环计算A的第i行与B的第j列的点积
3. 将结果写入C[i][j]

特点分析

优点：

• ✓ 代码简单，易于理解
• ✓ 结果完全正确
• ✓ 不需要特殊硬件

缺点：

• ✗ 完全串行，无法利用多核
• ✗ 性能较差

性能数据

在典型的现代CPU上（如Intel i7），性能表现如下：

可以看到，当矩阵规模翻倍时，计算时间增加了约8倍，这正好符合O(N³)的复杂度。

五、版本2：GPU Naive实现

现在，让我们将计算移到GPU上，利用并行计算的力量。

核心思想

GPU Naive版本的思路很简单：

• 每个thread计算C矩阵的一个元素
• 直接从Global Memory读取A和B的数据
• 并行执行，但尚未进行内存优化

代码实现

__global__ void matmul_naive(const float* A, const float* B, 
                             float* C, int M, int N, int K) {
    // 计算当前thread负责的元素位置
    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;
    }
}

Host端调用代码

void launch_matmul_naive(const float* A, const float* B, float* C,
                         int M, int N, int K) {
    // 定义block大小（通常为16×16或32×32）
    dim3 blockDim(16, 16);
    
    // 计算需要的block数量
    dim3 gridDim((N + blockDim.x - 1) / blockDim.x,
                 (M + blockDim.y - 1) / blockDim.y);
    
    // 启动kernel
    matmul_naive<<<gridDim, blockDim>>>(A, B, C, M, N, K);
}

内存访问分析

让我们仔细分析这个版本的内存访问模式：

计算C[row][col]需要：
┌─────────────────────────────────────┐
│ 读A[row][0...K-1]: K次Global Memory │
│ 读B[0...K-1][col]: K次Global Memory │
│ 写C[row][col]: 1次Global Memory     │
└─────────────────────────────────────┘

问题：
❌ 每次读取都访问Global Memory（延迟~400 cycles）
❌ A的每一行被不同的threads重复读取N次
❌ B的每一列被不同的threads重复读取M次
→ 存在大量重复的Global Memory访问！

虽然计算是并行的，但内存访问效率很低。这就是我们下一个版本要解决的问题。

性能数据对比

性能分析

从数据可以看出：

• GPU版本比CPU快了10-40倍，这主要来自并行计算的优势
• 矩阵越大，加速比越高（因为更好地利用了GPU的并行能力）
• 但这还远远不够！cuBLAS等高度优化的库能达到100倍以上的加速

为什么还有这么大的优化空间？因为我们还没有充分利用GPU的内存层级结构。接下来的Tiled版本将解决这个问题。

六、版本3：GPU Tiled优化实现 ⭐⭐⭐

这是本文的核心部分。通过使用Shared Memory和Tiling技术，我们将看到性能的巨大提升。

什么是Tiling？

Tiling（分块） 是一种经典的优化技术，核心思想是：将大矩阵分成小块（tiles/瓦片），每次只处理一个小块。

原始大矩阵（难以放入Shared Memory）:
┌─────────────────────────┐
│                         │
│    完整的1024×1024      │
│        矩阵             │
│                         │
└─────────────────────────┘

Tiling后（每个小块可以放入Shared Memory）:
┌───┬───┬───┬───┬───┬───┐
│T11│T12│T13│T14│T15│T16│  ← 每个Tile是16×16
├───┼───┼───┼───┼───┼───┤
│T21│T22│T23│T24│T25│T26│
├───┼───┼───┼───┼───┼───┤
│T31│T32│T33│T34│T35│T36│
└───┴───┴───┴───┴───┴───┘

策略：每次将一个Tile加载到Shared Memory中处理！

为什么Tiling能提升性能？

让我们通过一个具体的例子来理解数据复用：

Naive版本的问题：

计算C的一个16×16 block需要：
- 从A读取：16 rows × K elements = 16K 次Global Memory读取
- 从B读取：K rows × 16 elements = 16K 次Global Memory读取
- 总计：32K 次Global Memory访问

每个数据读一次，用一次，完全没有复用！

Tiled版本的优势：

将K维度分成 K/16 个tiles，每个tile处理：

步骤1：加载A的16×16 tile到Shared Memory（256次读取）
步骤2：加载B的16×16 tile到Shared Memory（256次读取）
步骤3：在Shared Memory中计算（16×16次乘加，无Global Memory访问）
步骤4：重复K/16次

总Global Memory读取：(K/16) × (256 + 256) = K × 32

关键：每个数据从Global Memory读一次，但在Shared Memory中被使用16次！
数据复用率：16倍
理论加速：接近16倍！

算法流程

Tiled矩阵乘法的完整流程如下：

对于C矩阵的每个block:
│
├─ 初始化：sum = 0
│
├─ for (tile_idx = 0; tile_idx < K/TILE_SIZE; tile_idx++)
│   │
│   ├─ 步骤1：协作加载A的当前tile到Shared Memory
│   │         (每个thread负责加载一个元素)
│   │
│   ├─ 步骤2：协作加载B的当前tile到Shared Memory
│   │         (每个thread负责加载一个元素)
│   │
│   ├─ 步骤3：__syncthreads() 
│   │         (确保所有数据都加载完成)
│   │
│   ├─ 步骤4：计算部分结果
│   │         (从Shared Memory读取，快速计算)
│   │
│   └─ 步骤5：__syncthreads()
│             (确保所有threads都用完数据，再加载新数据)
│
└─ 写入最终结果到Global Memory

完整代码实现

#define TILE_SIZE 16

__global__ void matmul_tiled(const float* A, const float* B, 
                             float* C, int M, int N, int K) {
    // ============================================
    // 步骤1：声明Shared Memory
    // ============================================
    // 两个tile，分别存储A和B的子矩阵
    __shared__ float tile_A[TILE_SIZE][TILE_SIZE];
    __shared__ float tile_B[TILE_SIZE][TILE_SIZE];
    
    // ============================================
    // 步骤2：计算thread的全局位置
    // ============================================
    int row = blockIdx.y * TILE_SIZE + threadIdx.y;
    int col = blockIdx.x * TILE_SIZE + threadIdx.x;
    
    // 累加器，存储最终结果
    float sum = 0.0f;
    
    // ============================================
    // 步骤3：遍历所有tiles
    // ============================================
    int num_tiles = (K + TILE_SIZE - 1) / TILE_SIZE;
    
    for (int t = 0; t < num_tiles; t++) {
        // ========================================
        // 步骤3.1：加载A的tile到Shared Memory
        // ========================================
        // 每个thread负责加载一个元素
        if (row < M && t * TILE_SIZE + threadIdx.x < K) {
            tile_A[threadIdx.y][threadIdx.x] = 
                A[row * K + t * TILE_SIZE + threadIdx.x];
        } else {
            // 边界外填充0（处理矩阵大小不是TILE_SIZE倍数的情况）
            tile_A[threadIdx.y][threadIdx.x] = 0.0f;
        }
        
        // ========================================
        // 步骤3.2：加载B的tile到Shared Memory
        // ========================================
        if (col < N && t * TILE_SIZE + threadIdx.y < K) {
            tile_B[threadIdx.y][threadIdx.x] = 
                B[(t * TILE_SIZE + threadIdx.y) * N + col];
        } else {
            // 边界外填充0
            tile_B[threadIdx.y][threadIdx.x] = 0.0f;
        }
        
        // ========================================
        // 步骤3.3：同步，确保tile加载完成
        // ========================================
        __syncthreads();
        
        // ========================================
        // 步骤3.4：计算部分结果
        // ========================================
        // 这里从Shared Memory读取，速度快！
        for (int k = 0; k < TILE_SIZE; k++) {
            sum += tile_A[threadIdx.y][k] * tile_B[k][threadIdx.x];
        }
        
        // ========================================
        // 步骤3.5：再次同步，确保计算完成
        // ========================================
        // 避免下一轮加载覆盖当前还在使用的数据
        __syncthreads();
    }
    
    // ============================================
    // 步骤4：写入最终结果
    // ============================================
    if (row < M && col < N) {
        C[row * N + col] = sum;
    }
}

Host端调用

void launch_matmul_tiled(const float* A, const float* B, float* C,
                         int M, int N, int K) {
    // Block大小必须与TILE_SIZE匹配
    dim3 blockDim(TILE_SIZE, TILE_SIZE);
    
    // 计算Grid大小
    dim3 gridDim((N + TILE_SIZE - 1) / TILE_SIZE,
                 (M + TILE_SIZE - 1) / TILE_SIZE);
    
    // 启动kernel
    matmul_tiled<<<gridDim, blockDim>>>(A, B, C, M, N, K);
}

关键技术点详解

1. 为什么需要两次__syncthreads()？

这是初学者最容易犯错的地方。让我们详细解释：

第一次同步（加载后）：

// Thread 0加载tile_A[0][0]
tile_A[0][0] = ...;

// 如果没有同步，Thread 1可能立即读取
// 但此时Thread 0可能还没写完！
__syncthreads();  // 必须等所有threads加载完

// 现在安全了，所有数据都ready
sum += tile_A[...][...] * tile_B[...][...];

第二次同步（计算后）：

// Thread 0还在使用tile_A
sum += tile_A[0][0] * ...;

// 如果没有同步，Thread 1可能开始加载新tile
// 覆盖了Thread 0还在用的数据！
__syncthreads();  // 必须等所有threads用完

// 现在安全了，可以加载新数据
tile_A[0][0] = new_value;

形象比喻：

• 第一次同步：等所有人把菜都端上桌，再开始吃
• 第二次同步：等所有人都吃完，再收拾桌子准备下一道菜

2. 边界检查的重要性

// 为什么需要边界检查？
if (row < M && t * TILE_SIZE + threadIdx.x < K) {
    tile_A[threadIdx.y][threadIdx.x] = A[...];
} else {
    tile_A[threadIdx.y][threadIdx.x] = 0.0f;  // 填充0
}

问题场景：
假设矩阵大小是1000×1000，TILE_SIZE=16

• 1000 / 16 = 62.5，需要63个tiles
• 最后一个tile只需要8个元素（1000 % 16 = 8）
• 剩余8个位置怎么办？填充0！

如果不填充0，这些位置是随机值，会导致计算结果错误。

3. 为什么选择TILE_SIZE=16？

这是一个权衡多个因素的结果：

Shared Memory容量限制：

每个block的Shared Memory：
tile_A: 16×16 × 4 bytes = 1024 bytes
tile_B: 16×16 × 4 bytes = 1024 bytes
总计：2048 bytes = 2 KB

GPU的Shared Memory：48 KB/block
2 KB << 48 KB，还有很大余量

太小的问题（如8×8）：

• 数据复用率只有8倍（而不是16倍）
• 优化效果不明显

太大的问题（如32×32）：

• Shared Memory使用：32×32×2×4 = 8 KB
• 限制了每个SM的active blocks数量
• 可能降低occupancy

经验值：

• 16×16：通用，适合大多数情况
• 32×32：适合大矩阵，需要更多Shared Memory
• 8×8：适合Shared Memory紧张或小矩阵的情况

七、性能对比与深度分析

终于到了揭晓答案的时刻！让我们看看三个版本的性能对比。

完整性能测试结果

测试环境：NVIDIA RTX 4060 (8GB), CUDA 12.8

性能提升可视化

性能对比 (1024×1024矩阵):

CPU:    ████████████████████████████████████████ 3600ms
        
Naive:  ████████ 180ms (20x faster than CPU)
        
Tiled:  ▌ 15ms (240x faster than CPU, 12x faster than Naive)

深度分析：为什么Tiled这么快？

让我们从多个角度分析性能提升的原因。

1. 内存访问次数大幅减少

定量分析（以1024×1024矩阵为例）：

Naive版本：
- 每个元素：2K 次Global Memory读取
- 总计：1024² × 2 × 1024 = 2,147,483,648 次读取
- 约21亿次Global Memory访问

Tiled版本（TILE_SIZE=16）：
- 每个tile：256 + 256 = 512 次Global Memory读取
- Tiles数量：(1024/16)² × (1024/16) = 262,144 个tiles
- 总计：262,144 × 512 = 134,217,728 次读取
- 约1.3亿次Global Memory访问

减少比例：2,147,483,648 / 134,217,728 ≈ 16倍！

2. 数据复用率提升

数据复用分析：

Naive：每个数据读一次，用一次
复用率 = 1

Tiled：每个数据从Global读一次，在Shared Memory中被用16次
复用率 = 16

提升：16倍！

3. 内存带宽利用率

理论分析：

Naive版本：
- 访问模式：随机、分散
- Cache命中率：低
- 带宽利用：~20-30%

Tiled版本：
- 访问模式：连续、块状
- Cache命中率：高
- 带宽利用：~60-70%

实际测量（使用nvprof）：
- Naive：内存带宽 150 GB/s
- Tiled：内存带宽 350 GB/s

4. 计算强度（Arithmetic Intensity）

这是一个衡量算法效率的重要指标：

AI = 浮点运算次数 / 内存访问字节数

Naive版本：
- FLOPs = 2K (每个元素K次乘法+K次加法)
- Memory = 2K × 4 bytes = 8K bytes
- AI = 2K / 8K = 0.25 FLOPs/Byte

Tiled版本：
- FLOPs = 2K (相同)
- Memory = (2K/16) × 4 bytes = K/2 bytes
- AI = 2K / (K/2) = 4 FLOPs/Byte

AI提升：4 / 0.25 = 16倍！

计算强度越高，说明每读取一个字节的数据，进行的计算越多，效率越高。

与cuBLAS的对比

cuBLAS是NVIDIA官方的高度优化的线性代数库。让我们看看我们的实现与它的差距：

cuBLAS比我们的Tiled版本还快12倍左右！这说明还有巨大的优化空间。cuBLAS使用了：

• 更复杂的tiling策略（3D tiling）
• Tensor Core硬件加速
• 向量化内存访问（float4）
• Bank conflict优化
• Warp-level primitives
• 汇编级优化

这也是为什么实际项目中应该使用成熟的库，而不是自己实现。

八、实战经验与进阶技巧

通过实现这三个版本，我们积累了一些宝贵的经验。让我分享一些实用技巧。

选择合适的TILE_SIZE

选择TILE_SIZE需要权衡多个因素：

权衡因素表

推荐策略

// 根据GPU架构选择
#if __CUDA_ARCH__ >= 750  // Turing及更新
    #define TILE_SIZE 32
#elif __CUDA_ARCH__ >= 600  // Pascal及更新
    #define TILE_SIZE 16
#else
    #define TILE_SIZE 8
#endif

避免Bank Conflict

Shared Memory被组织成32个banks。如果同一warp的多个threads访问同一bank的不同地址，会发生bank conflict，导致访问串行化。

问题示例

__shared__ float tile[16][16];

// 可能存在bank conflict
for (int i = 0; i < 16; i++) {
    sum += tile[i][threadIdx.x];  // 所有threads访问同一列
}

解决方案：Padding

// 添加padding避免bank conflict
__shared__ float tile[16][17];  // 注意：17 = 16 + 1

// 现在访问被打散了，减少conflict
for (int i = 0; i < 16; i++) {
    sum += tile[i][threadIdx.x];
}

完整的性能优化checklist

在实际项目中，你可以按照这个checklist逐步优化：

□ Level 1：基础并行化
  ├─ □ 使用GPU替代CPU
  ├─ □ 合理设置block和grid大小
  └─ □ 确保正确性

□ Level 2：内存优化
  ├─ □ 使用Shared Memory
  ├─ □ 实现Tiling
  ├─ □ 优化内存访问模式
  └─ □ 添加边界检查

□ Level 3：细节优化
  ├─ □ 选择最优TILE_SIZE
  ├─ □ 避免bank conflict
  ├─ □ 使用向量化加载(float4)
  └─ □ 优化寄存器使用

□ Level 4：高级优化
  ├─ □ 使用Tensor Cores
  ├─ □ 实现double buffering
  ├─ □ Warp-level primitives
  └─ □ 汇编级优化

□ Level 5：验证与测试
  ├─ □ 性能测试（多次运行取平均）
  ├─ □ 正确性验证
  ├─ □ 使用nvprof分析
  └─ □ 与cuBLAS对比

完整代码（matmul.cu）

#include <iostream>
#include <cuda_runtime.h>
#include <cmath>
#include <chrono>

#define TILE_SIZE 16

// ============================================
// 版本1：CPU实现
// ============================================
void matmul_cpu(const float* A, const float* B, float* C, 
                int M, int N, int K) {
    for (int i = 0; i < M; i++) {
        for (int j = 0; j < N; j++) {
            float sum = 0.0f;
            for (int k = 0; k < K; k++) {
                sum += A[i * K + k] * B[k * N + j];
            }
            C[i * N + j] = sum;
        }
    }
}

// ============================================
// 版本2：GPU Naive实现
// ============================================
__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;
    }
}

// ============================================
// 版本3：GPU Tiled实现
// ============================================
__global__ void matmul_tiled(const float* A, const float* B, 
                             float* C, int M, int N, int K) {
    __shared__ float tile_A[TILE_SIZE][TILE_SIZE];
    __shared__ float tile_B[TILE_SIZE][TILE_SIZE];
    
    int row = blockIdx.y * TILE_SIZE + threadIdx.y;
    int col = blockIdx.x * TILE_SIZE + threadIdx.x;
    
    float sum = 0.0f;
    int num_tiles = (K + TILE_SIZE - 1) / TILE_SIZE;
    
    for (int t = 0; t < num_tiles; t++) {
        if (row < M && t * TILE_SIZE + threadIdx.x < K) {
            tile_A[threadIdx.y][threadIdx.x] = 
                A[row * K + t * TILE_SIZE + threadIdx.x];
        } else {
            tile_A[threadIdx.y][threadIdx.x] = 0.0f;
        }
        
        if (col < N && t * TILE_SIZE + threadIdx.y < K) {
            tile_B[threadIdx.y][threadIdx.x] = 
                B[(t * TILE_SIZE + threadIdx.y) * N + col];
        } else {
            tile_B[threadIdx.y][threadIdx.x] = 0.0f;
        }
        
        __syncthreads();
        
        for (int k = 0; k < TILE_SIZE; k++) {
            sum += tile_A[threadIdx.y][k] * tile_B[k][threadIdx.x];
        }
        
        __syncthreads();
    }
    
    if (row < M && col < N) {
        C[row * N + col] = sum;
    }
}

// ============================================
// 辅助函数
// ============================================
void init_matrix(float* mat, int size) {
    for (int i = 0; i < size; i++) {
        mat[i] = static_cast<float>(rand()) / RAND_MAX;
    }
}

bool verify(const float* C1, const float* C2, int size, float eps = 1e-3) {
    for (int i = 0; i < size; i++) {
        if (fabs(C1[i] - C2[i]) > eps) {
            std::cout << "Mismatch at " << i << ": " 
                      << C1[i] << " vs " << C2[i] << std::endl;
            return false;
        }
    }
    return true;
}

// ============================================
// 主函数
// ============================================
int main() {
    const int M = 1024;
    const int N = 1024;
    const int K = 1024;
    
    // 分配host内存
    float *h_A, *h_B, *h_C_cpu, *h_C_naive, *h_C_tiled;
    h_A = new float[M * K];
    h_B = new float[K * N];
    h_C_cpu = new float[M * N];
    h_C_naive = new float[M * N];
    h_C_tiled = new float[M * N];
    
    // 初始化
    init_matrix(h_A, M * K);
    init_matrix(h_B, K * N);
    
    // 分配device内存
    float *d_A, *d_B, *d_C;
    cudaMalloc(&d_A, M * K * sizeof(float));
    cudaMalloc(&d_B, K * N * sizeof(float));
    cudaMalloc(&d_C, M * N * sizeof(float));
    
    cudaMemcpy(d_A, h_A, M * K * sizeof(float), cudaMemcpyHostToDevice);
    cudaMemcpy(d_B, h_B, K * N * sizeof(float), cudaMemcpyHostToDevice);
    
    // 测试CPU版本
    std::cout << "Testing CPU version..." << std::endl;
    auto start = std::chrono::high_resolution_clock::now();
    matmul_cpu(h_A, h_B, h_C_cpu, M, N, K);
    auto end = std::chrono::high_resolution_clock::now();
    auto cpu_time = std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count();
    
    // 测试Naive版本
    std::cout << "Testing Naive GPU version..." << std::endl;
    dim3 blockDim(16, 16);
    dim3 gridDim((N + 15) / 16, (M + 15) / 16);
    
    cudaEvent_t start_event, stop_event;
    cudaEventCreate(&start_event);
    cudaEventCreate(&stop_event);
    
    cudaEventRecord(start_event);
    matmul_naive<<<gridDim, blockDim>>>(d_A, d_B, d_C, M, N, K);
    cudaEventRecord(stop_event);
    cudaEventSynchronize(stop_event);
    
    float naive_time;
    cudaEventElapsedTime(&naive_time, start_event, stop_event);
    cudaMemcpy(h_C_naive, d_C, M * N * sizeof(float), cudaMemcpyDeviceToHost);
    
    // 测试Tiled版本
    std::cout << "Testing Tiled GPU version..." << std::endl;
    dim3 gridDim2((N + TILE_SIZE - 1) / TILE_SIZE, 
                  (M + TILE_SIZE - 1) / TILE_SIZE);
    
    cudaEventRecord(start_event);
    matmul_tiled<<<gridDim2, blockDim>>>(d_A, d_B, d_C, M, N, K);
    cudaEventRecord(stop_event);
    cudaEventSynchronize(stop_event);
    
    float tiled_time;
    cudaEventElapsedTime(&tiled_time, start_event, stop_event);
    cudaMemcpy(h_C_tiled, d_C, M * N * sizeof(float), cudaMemcpyDeviceToHost);
    
    // 验证正确性
    std::cout << "\n=== Verification ===" << std::endl;
    std::cout << "Naive vs CPU: " 
              << (verify(h_C_cpu, h_C_naive, M * N) ? "PASS" : "FAIL") 
              << std::endl;
    std::cout << "Tiled vs CPU: " 
              << (verify(h_C_cpu, h_C_tiled, M * N) ? "PASS" : "FAIL") 
              << std::endl;
    
    // 性能报告
    std::cout << "\n=== Performance ===" << std::endl;
    std::cout << "Matrix size: " << M << "×" << N << std::endl;
    std::cout << "CPU time: " << cpu_time << " ms" << std::endl;
    std::cout << "Naive GPU time: " << naive_time << " ms ("
              << cpu_time / naive_time << "x faster)" << std::endl;
    std::cout << "Tiled GPU time: " << tiled_time << " ms ("
              << cpu_time / tiled_time << "x faster than CPU, "
              << naive_time / tiled_time << "x faster than Naive)" << std::endl;
    
    // 清理
    delete[] h_A;
    delete[] h_B;
    delete[] h_C_cpu;
    delete[] h_C_naive;
    delete[] h_C_tiled;
    cudaFree(d_A);
    cudaFree(d_B);
    cudaFree(d_C);
    cudaEventDestroy(start_event);
    cudaEventDestroy(stop_event);
    
    return 0;
}

预期输出

Testing CPU version...
Testing Naive GPU version...
Testing Tiled GPU version...

=== Verification ===
Naive vs CPU: PASS
Tiled vs CPU: PASS

=== Performance ===
Matrix size: 1024×1024
CPU time: 3600 ms
Naive GPU time: 180 ms (20.0x faster)
Tiled GPU time: 15 ms (240.0x faster than CPU, 12.0x faster than Naive)

十、总结与展望

经过这篇文章的学习，我们完成了从CPU到GPU，从Naive到Tiled的完整优化旅程。

核心收获总结

理论知识

我们掌握了以下核心概念：

• ✓ GPU内存层级：理解了Registers、Shared Memory和Global Memory的特性
• ✓ Shared Memory优化：学会了如何使用Shared Memory加速计算
• ✓ Tiling思想：掌握了分块处理的经典优化技术
• ✓ 数据复用：理解了为什么数据复用是优化的关键

实战能力

我们获得了以下技能：

• ✓ 三种实现方式：CPU、GPU Naive、GPU Tiled
• ✓ 性能测试方法：使用cudaEvent进行精确计时
• ✓ 正确性验证：对比不同版本的计算结果
• ✓ 性能分析：理解性能提升的原因

关键Insights

核心insight #1: 内存访问是性能瓶颈
→ Global Memory延迟是Shared Memory的200倍

核心insight #2: 数据复用是优化的本质
→ Tiling将复用率从1提升到16

核心insight #3: 同步是正确性的保证
→ 两次__syncthreads()缺一不可

核心insight #4: 理论分析与实测结合
→ 计算强度、带宽利用率等指标指导优化

Shared Memory优化的通用模式

通过矩阵乘法的案例，我们提炼出Shared Memory优化的通用模式：

Step 1: 分析数据访问模式
        └─ 识别哪些数据会被多次访问

Step 2: 设计Tiling策略
        └─ 确定合适的tile大小

Step 3: 协作加载数据
        └─ 每个thread加载一部分数据到Shared Memory

Step 4: 同步等待
        └─ __syncthreads()确保数据ready

Step 5: 快速计算
        └─ 从Shared Memory读取，高速计算

Step 6: 再次同步
        └─ __syncthreads()确保数据用完

Step 7: 写回结果
        └─ 将最终结果写入Global Memory

何时使用Shared Memory？

适用场景：

• ✓ 数据会被Block内多个threads重复读取
• ✓ Threads之间需要协作和通信
• ✓ 内存访问是性能瓶颈（而非计算）
• ✓ 数据大小适合放入Shared Memory

不适用场景：

• ✗ 每个数据只用一次
• ✗ Threads之间完全独立
• ✗ 计算本身很慢（计算密集型）
• ✗ 数据量太大，无法放入Shared Memory

性能优化的层次

我们的优化之旅可以总结为：

Level 1: CPU串行 (基准)
         └─ 简单，但慢

Level 2: GPU并行 (20x)
         └─ 利用并行计算

Level 3: Shared Memory (240x)
         └─ 优化内存访问

Level 4: 高级优化 (cuBLAS水平，1000x+)
         ├─ Double buffering
         ├─ Vectorized load (float4)
         ├─ Bank conflict优化
         ├─ Warp-level primitives
         └─ Tensor Cores

每一层都带来数量级的提升！