第4章 GPU内存架构 - 显存的层次与优化 (GPU Memory Architecture - Hierarchy and Optimization)

一、导读

4.1 本章学习目标

通过本章的学习，读者将能够：

全面理解GPU内存层次的组成：寄存器、共享内存、L1/L2缓存、全局显存、常量内存、纹理内存

掌握每种内存类型的特性（速度、容量、作用域、生命周期）及其适用场景

深入理解全局显存的合并访问（coalesced access）机制，能够优化内存访问模式

掌握共享内存的编程技巧，包括数据复用、线程间通信和bank conflict避免

理解L1/L2缓存的工作原理，以及如何利用缓存提高内存访问效率

掌握常量内存和纹理内存的特殊用途及优化机制

学会使用Nsight Compute分析内存带宽利用率和缓存命中率

能够根据应用特点设计最优的数据布局和内存访问策略

4.2 为什么本章如此重要

在GPU编程中，内存性能往往是决定整体性能的关键瓶颈。与现代CPU类似，GPU也采用了多级内存层次结构来平衡速度、容量和成本。然而，GPU的内存模型与CPU有显著差异：GPU更强调程序员对内存的显式管理，而非完全依赖硬件缓存。

理解GPU内存架构对于编写高性能CUDA代码至关重要。许多初学者编写的kernel计算逻辑正确，但性能远低于预期，主要原因就是内存访问模式不佳。例如，非合并的全局内存访问可能导致带宽利用率下降一个数量级；共享内存的bank conflict可能使并行访问退化为串行访问；过度依赖全局内存而忽视数据复用会导致大量的冗余数据传输。

GPU的内存带宽虽然远高于CPU主存，但仍然有限。Ampere架构的HBM2e显存带宽约为1.6 TB/s，看似很高，但当数千个线程同时访问内存时，每个线程能分到的带宽非常有限。因此，优化内存访问是GPU性能调优的核心任务之一。

此外，不同的内存类型有不同的编程模型和使用约束。寄存器最快但容量有限；共享内存速度快且可程序控制，但需要显式管理同步；全局显存容量大但速度慢；常量内存和纹理内存针对特定访问模式进行了优化。只有深入理解每种内存的特点，才能在实际应用中做出正确的选择。

4.3 核心问题引导

GPU的内存层次结构是怎样的？每种内存的速度、容量和作用域如何？

什么是合并访问（coalesced access）？如何实现合并访问？

共享内存与L1缓存有什么区别？何时应该使用共享内存而非依赖缓存？

什么是bank conflict？如何检测和避免bank conflict？

常量内存的广播机制是如何工作的？它适用于哪些场景？

纹理内存的缓存优势是什么？除了图形应用，它在通用计算中有哪些用途？

L2缓存在多SM系统中扮演什么角色？如何利用L2缓存提高数据共享效率？

统一内存（Unified Memory）与传统显存管理有何不同？它的优势和局限是什么？

二、核心概念详解

4.1 GPU内存层次概览

GPU的内存层次结构从快到慢、从小到大依次为：

内存类型	位置	访问延迟	带宽	容量	作用域	管理方式
寄存器	SM内部	~1周期	极高	64K/SM	线程私有	编译器分配
共享内存	SM内部	~1-2周期	极高	100KB/SM	Block内共享	程序员管理
L1缓存	SM内部	~20周期	高	可变	线程/Block	硬件管理
L2缓存	芯片级	~100周期	中	MB级别	全局	硬件管理
全局显存	板载	~400-800周期	高	GB级别	全局	程序员分配
常量内存	专用缓存	~1周期(命中)	低	64KB	只读	程序员管理
纹理内存	专用缓存	~1周期(命中)	低	可变	只读	程序员管理

寄存器（Registers）：

寄存器是GPU最快的存储单元，位于SM内部的寄存器文件中。每个线程有自己的寄存器空间，用于存储局部变量和中间计算结果。寄存器的访问延迟仅为1个时钟周期，带宽极高，因为每个线程可以直接访问自己的寄存器而无需仲裁。

寄存器的主要限制是容量有限。Ampere SM的寄存器文件包含65536个32位寄存器，如果每个线程使用较多寄存器，则会限制SM可以同时驻留的线程数量，影响Occupancy。编译器负责将源代码中的变量映射到物理寄存器，当寄存器不足时，会将部分变量溢出到本地内存（实际上是全局内存的一部分），导致性能急剧下降。

共享内存（Shared Memory）：

共享内存是GPU特有的高速片上存储器，位于SM内部。它的速度与寄存器相当，但可以被同一Block内的所有线程共享。共享内存的主要用途包括：

线程间通信：Block内的线程可以通过共享内存交换数据

数据复用：将频繁访问的数据加载到共享内存，减少全局内存访问

软件管理的缓存：程序员可以精确控制哪些数据放入共享内存，比硬件缓存更可预测

共享内存由程序员显式管理，需要使用__shared__关键字声明，并通过__syncthreads()确保数据一致性。共享内存被划分为多个bank，需要注意bank conflict问题。

L1缓存：

L1缓存是位于SM内部的第一级缓存，主要用于加速全局内存和局部内存的访问。在较新的GPU架构中，L1缓存与共享内存共享物理存储空间，可以动态配置两者的比例。L1缓存由硬件自动管理，程序员无法直接控制数据的缓存行为。

L1缓存的优势是透明性：程序员无需修改代码即可享受缓存带来的加速。但缺点是可控性差：对于不规则的访问模式，缓存命中率可能很低。

L2缓存：

L2缓存是芯片级的二级缓存，所有SM共享。L2缓存的容量通常为几MB到几十MB，远大于L1缓存。L2缓存的作用包括：

跨SM数据共享：不同SM上的线程可以通过L2缓存共享数据

减少全局内存访问：L2缓存可以捕获重复的全局内存访问

写缓冲：写操作可以先写入L2缓存，再异步写回全局内存

L2缓存同样由硬件管理，对所有内存访问透明。

全局显存（Global Memory）：

全局显存是GPU上容量最大的存储器，通常为数GB到数十GB。全局显存使用GDDR6或HBM技术，带宽远高于CPU主存，但访问延迟也较高（数百个时钟周期）。全局内存是所有线程都可以访问的唯一存储器，也是主机与设备之间数据传输的目标。

全局内存的性能高度依赖于访问模式。合并访问可以实现接近理论峰值的带宽，而非合并访问可能导致带宽利用率下降一个数量级。

常量内存（Constant Memory）：

常量内存是专门用于存储只读常量的存储器，容量为64KB。常量内存通过专用的常量缓存进行访问，当warp内所有线程读取相同的常量地址时，硬件可以将单次内存事务广播给所有线程，实现极高的效率。

常量内存适用于存储kernel参数、查找表等不变的数据。使用__constant__关键字声明常量内存变量。

纹理内存（Texture Memory）：

纹理内存最初为图形渲染设计，用于存储纹理数据。纹理内存通过专用的纹理缓存进行访问，支持硬件级别的插值、地址模式转换等操作。在通用计算中，纹理内存可以用于具有空间局部性的数据访问模式，如图像处理、卷积计算等。

纹理内存的访问通过特殊的纹理获取指令完成，可以利用二维空间局部性提高缓存命中率。

4.2 全局内存与合并访问

全局内存是GPU上最主要的存储器，但其性能高度依赖于访问模式。理解合并访问机制是优化全局内存访问的关键。

合并访问的定义：

合并访问是指warp内的线程访问全局内存时，它们的地址落在同一个或少数几个连续的内存段（segment）内，使得硬件可以用最少的内存事务完成所有线程的访问请求。

在NVIDIA GPU中，全局内存访问以32字节、64字节或128字节的段为单位进行。理想情况下，warp内的32个线程应该访问连续的128字节数据（每个线程4字节），这样只需要一次128字节的内存事务即可完成。

合并访问的条件：

对于32位数据类型（如float、int），合并访问的条件是：

warp内的线程访问连续的内存地址

起始地址是适当对齐的（通常是128字节对齐）

具体来说，如果thread 0访问地址A，thread 1访问地址A+4，thread 2访问地址A+8，...，thread 31访问地址A+124，则这是一个完美的合并访问，只需要一次128字节的内存事务。

非合并访问的例子：

cpp

// Non-coalesced: threads access strided addresses
__global__ void badAccess(float* data, int N) {
    int idx = blockIdx.x * blockDim.x + threadIdx.x;
    if (idx < N) {
        // Thread 0 accesses data[0], thread 1 accesses data[32], etc.
        data[idx * 32] = compute(idx);
    }
}

// Coalesced: threads access consecutive addresses
__global__ void goodAccess(float* data, int N) {
    int idx = blockIdx.x * blockDim.x + threadIdx.x;
    if (idx < N) {
        // Thread 0 accesses data[0], thread 1 accesses data[1], etc.
        data[idx] = compute(idx);
    }
}

在badAccess中，相邻线程访问的地址相差128字节（32个float），导致每个线程都需要独立的内存事务，带宽利用率极低。在goodAccess中，相邻线程访问连续地址，实现了完美合并。

转置访问的问题：

矩阵按行存储在内存中（行优先），如果按列访问矩阵，会导致非合并访问：

cpp

// Column-major access on row-major stored matrix: non-coalesced
__global__ void columnAccess(float* matrix, int width, int height) {
    int col = blockIdx.x * blockDim.x + threadIdx.x;
    int row = blockIdx.y * blockDim.y + threadIdx.y;
    
    if (col < width && row < height) {
        // Accessing by column: stride = width
        float val = matrix[row * width + col];
        // If adjacent threads have different row but same col, 
        // they access addresses far apart
    }
}

解决这个问题的方法是使用共享内存进行数据重排，或者改变算法的访问模式。

合并访问的量化指标：

Nsight Compute工具可以报告以下内存效率指标：

Memory Throughput：实际达到的内存带宽
Compute Throughput：实际达到的计算吞吐量
L1/L2 Cache Hit Rate：缓存命中率
DRAM Utilization：显存带宽利用率

理想的合并访问应该达到80%以上的DRAM利用率。

4.3 共享内存编程技巧

共享内存是GPU编程中最强大的优化工具之一，正确使用共享内存可以显著提升性能。

数据复用模式：

当多个线程需要访问相同的全局内存数据时，可以先将数据加载到共享内存，然后从共享内存读取：

cpp

__global__ void reusePattern(const float* input, float* output, int N) {
    extern __shared__ float shared[];
    
    int tid = threadIdx.x;
    int idx = blockIdx.x * blockDim.x + threadIdx.x;
    
    // Load data into shared memory
    if (idx < N) {
        shared[tid] = input[idx];
    }
    __syncthreads();
    
    // Multiple threads can reuse the same shared data
    float result = 0.0f;
    for (int i = 0; i < windowSize; i++) {
        int neighborIdx = (tid + i) % blockDim.x;
        result += shared[neighborIdx];
    }
    
    if (idx < N) {
        output[idx] = result;
    }
}

这种模式在卷积、滤波等应用中非常常见。

矩阵乘法中的共享内存优化：

经典的矩阵乘法优化示例：

cpp

#define TILE_SIZE 16

__global__ void matMulShared(const float* A, const float* B, float* C,
                              int M, int N, int K) {
    __shared__ float As[TILE_SIZE][TILE_SIZE];
    __shared__ float Bs[TILE_SIZE][TILE_SIZE];
    
    int tx = threadIdx.x;
    int ty = threadIdx.y;
    int row = blockIdx.y * TILE_SIZE + ty;
    int col = blockIdx.x * TILE_SIZE + tx;
    
    float sum = 0.0f;
    
    // Loop over tiles
    for (int t = 0; t < (K + TILE_SIZE - 1) / TILE_SIZE; t++) {
        // Load tile of A into shared memory
        if (row < M && t * TILE_SIZE + tx < K) {
            As[ty][tx] = A[row * K + t * TILE_SIZE + tx];
        } else {
            As[ty][tx] = 0.0f;
        }
        
        // Load tile of B into shared memory
        if (t * TILE_SIZE + ty < K && col < N) {
            Bs[ty][tx] = B[(t * TILE_SIZE + ty) * N + col];
        } else {
            Bs[ty][tx] = 0.0f;
        }
        
        __syncthreads();
        
        // Compute partial result for this tile
        for (int k = 0; k < TILE_SIZE; k++) {
            sum += As[ty][k] * Bs[k][tx];
        }
        
        __syncthreads();
    }
    
    if (row < M && col < N) {
        C[row * N + col] = sum;
    }
}

这个实现将全局内存访问减少了TILE_SIZE倍，因为每个元素从全局内存加载一次，然后被TILE_SIZE个线程复用。

Bank Conflict避免：

共享内存被划分为32个bank，每个bank每个周期可以服务一个32位字的访问。如果多个线程访问同一个bank的不同地址，就会发生bank conflict，访问需要串行化。

避免bank conflict的常用技巧：

填充（Padding）：在数组第二维加1

cpp

// Causes bank conflict when tx accesses column tx
__shared__ float shared[32][32];
shared[ty][tx] = value;  // tx=0 and tx=32 access same bank

// No bank conflict with padding
__shared__ float shared[32][33];
shared[ty][tx] = value;  // Different banks due to padding

广播机制：如果warp内所有线程访问相同地址，不会发生bank conflict，硬件会广播该值。

数据重排：重新组织数据结构，使得相邻线程访问相邻地址。

4.4 常量内存与纹理内存

常量内存的广播机制：

常量内存通过专用的常量缓存访问，容量为64KB。当warp内所有线程读取相同的常量地址时，硬件可以将单次内存读取广播给所有32个线程，只需1个时钟周期。这使得常量内存非常适合存储kernel参数、系数表等不变的数据。

cpp

__constant__ float coefficients[10];

__global__ void useConstants(float* data, int N) {
    int idx = blockIdx.x * blockDim.x + threadIdx.x;
    if (idx < N) {
        float result = 0.0f;
        for (int i = 0; i < 10; i++) {
            result += data[idx + i] * coefficients[i];
        }
        data[idx] = result;
    }
}

在这个例子中，所有线程读取相同的coefficients[i]，常量内存的广播机制使得这些读取非常高效。

但如果warp内线程读取不同的常量地址，常量缓存的效率会下降，可能需要多次内存事务。

纹理内存的空间局部性优化：

纹理内存通过专用的纹理缓存访问，针对二维空间局部性进行了优化。纹理缓存使用空间块（spatial blocking）策略，当访问某个地址时，会预取周围的地址到缓存中。这对于图像处理、卷积等具有空间局部性的应用非常有利。

在CUDA中，可以使用CUDA Array或绑定纹理引用到线性内存来使用纹理内存：

cpp

texture<float, 2, cudaReadModeElementType> texRef;

__global__ void textureKernel(float* output, int width, int height) {
    int x = blockIdx.x * blockDim.x + threadIdx.x;
    int y = blockIdx.y * blockDim.y + threadIdx.y;
    
    if (x < width && y < height) {
        // Texture fetch with automatic caching
        float val = tex2D(texRef, x, y);
        output[y * width + x] = val * 2.0f;
    }
}

纹理内存还支持硬件级别的线性插值和归一化坐标，在图像缩放等应用中非常方便。

4.5 L1/L2缓存的行为与优化

L1缓存：

在Ampere架构中，L1缓存与共享内存共享192KB的物理存储空间。默认配置下，L1缓存为80KB，共享内存为112KB，但可以通过CUDA API调整比例。

L1缓存主要缓存全局内存和本地内存的访问。对于具有时间局部性的访问模式（即短时间内重复访问相同地址），L1缓存可以显著提高性能。但对于一次性访问或随机访问，L1缓存的命中率可能很低。

L2缓存：

L2缓存是所有SM共享的二级缓存，Ampere A100的L2缓存大小为40MB。L2缓存的作用包括：

跨SM数据共享：当不同SM上的线程访问相同的全局内存地址时，L2缓存可以提供数据，避免重复访问全局内存

写合并：多个写操作可以在L2缓存中合并，减少全局内存写入次数

原子操作加速：L2缓存可以加速全局内存上的原子操作

L2缓存的大小对性能有显著影响。较大的L2缓存可以捕获更多的工作集，减少全局内存访问。

缓存友好编程：

虽然缓存由硬件管理，但程序员可以通过以下方式提高缓存效率：

数据局部性：尽量访问相邻的数据，利用空间局部性

时间局部性：重复使用的数据尽量保持在缓存中

减少工作集：如果数据可以放入L2缓存，性能会显著提升

避免缓存污染：不要访问大量只用一次的数据，这会挤出有用的缓存行

4.6 统一内存（Unified Memory）

统一内存是CUDA 6.0引入的功能，允许主机和设备共享同一虚拟地址空间。程序员可以分配统一内存，然后在主机和设备代码中直接使用指针访问，无需显式的数据传输。

工作原理：

统一内存基于页面迁移机制。数据最初可能在主机内存或设备内存中，当另一方访问时，页面会自动迁移到访问者的物理内存中。这种按需迁移简化了编程模型，但可能带来性能开销。

优势：

编程简化：无需手动管理cudaMemcpy

超出显存容量的数据处理：可以访问超过GPU显存容量的数据，页面会按需迁移

稀疏或不规则数据访问：对于难以预先确定访问模式的应用，统一内存更加灵活

局限：

页面迁移开销：频繁的跨处理器访问会导致大量的页面迁移，性能下降

缺页异常处理：首次访问远程数据时会触发缺页异常，需要操作系统介入

不可预测的性能：由于页面迁移的动态性，性能可能不稳定

优化建议：

使用cudaMemAdvise提示数据的使用模式
使用cudaPrefetchAsync预先将数据迁移到目标处理器
对于性能敏感的代码，仍然建议使用显式的cudaMalloc和cudaMemcpy

三、程序员的视角

3.1 在实际项目中的应用

深度学习框架：cuDNN等库大量使用共享内存优化卷积和归一化操作。卷积核权重通常存储在常量内存或纹理内存中，利用广播和缓存优势。激活值和梯度在全局内存中，但通过共享内存暂存实现高效的数据复用。

科学计算：在分子动力学模拟中，粒子位置存储在全局内存，但邻近列表（neighbor list）可以加载到共享内存中，减少重复的全局内存访问。力计算涉及大量随机访问，L2缓存的大小对性能影响显著。

数据库加速：GPU加速的数据库连接操作使用共享内存构建哈希表，线程协作查询。常量内存用于存储查询参数和配置信息。

实时渲染：光线追踪中，BVH树节点通过纹理内存访问，利用空间局部性提高缓存命中率。材质参数存储在常量内存中，供所有着色器线程读取。

3.2 常见陷阱与误区

忽视合并访问：这是最常见的性能问题。程序员应该始终检查全局内存访问模式，确保相邻线程访问相邻地址。

过度使用共享内存：共享内存容量有限，过度使用会限制Occupancy。应该在数据复用收益和Occupancy损失之间权衡。

Bank Conflict忽视：共享内存的bank conflict可能导致严重的性能下降，但往往被忽视。使用Nsight Compute检测bank conflict。

常量内存滥用：当warp内线程访问不同的常量地址时，常量内存的优势消失。确保常量内存用于真正的常量数据，且访问模式一致。

统一内存的性能陷阱：统一内存简化了编程，但页面迁移开销可能很大。对于性能敏感的代码，建议使用显式内存管理。

本地内存溢出：当寄存器不足时，编译器会将变量溢出到本地内存（全局内存的一部分），导致性能急剧下降。监控寄存器使用量，必要时限制最大寄存器数。

3.3 与其他层级的联系

与线程层次的联系：共享内存的作用域是Block，理解Block的边界对于设计有效的数据复用策略至关重要。寄存器是线程私有的，适合存储临时变量。

与执行模型的联系：合并访问的效率依赖于warp的执行特性。理解warp如何发起内存请求有助于优化访问模式。

与性能分析的联系：Nsight Compute提供了详细的内存性能指标，包括带宽利用率、缓存命中率、bank conflict等。这些指标指导内存优化方向。

四、代码实战

4.1 合并访问演示与对比

cpp

#include <cuda_runtime.h>
#include <stdio.h>

// Non-coalesced access pattern
__global__ void nonCoalesced(float* output, const float* input, int N) {
    int idx = blockIdx.x * blockDim.x + threadIdx.x;
    if (idx < N) {
        // Strided access: thread i accesses element i * 32
        output[idx] = input[idx * 32 % N];
    }
}

// Coalesced access pattern
__global__ void coalesced(float* output, const float* input, int N) {
    int idx = blockIdx.x * blockDim.x + threadIdx.x;
    if (idx < N) {
        // Sequential access: thread i accesses element i
        output[idx] = input[idx];
    }
}

// Transposed access (common pitfall)
__global__ void transposedAccess(float* output, const float* input, 
                                  int width, int height) {
    int x = blockIdx.x * blockDim.x + threadIdx.x;
    int y = blockIdx.y * blockDim.y + threadIdx.y;
    
    if (x < width && y < height) {
        // Reading by column from row-major stored matrix
        output[y * width + x] = input[x * height + y];
    }
}

// Benchmark function
void benchmarkMemoryAccess(int N) {
    size_t size = N * sizeof(float);
    
    float *h_input = (float*)malloc(size);
    float *h_output = (float*)malloc(size);
    
    for (int i = 0; i < N; i++) {
        h_input[i] = (float)i;
    }
    
    float *d_input, *d_output;
    cudaMalloc((void**)&d_input, size);
    cudaMalloc((void**)&d_output, size);
    
    cudaMemcpy(d_input, h_input, size, cudaMemcpyHostToDevice);
    
    int blockSize = 256;
    int gridSize = (N + blockSize - 1) / blockSize;
    
    cudaEvent_t start, stop;
    cudaEventCreate(&start);
    cudaEventCreate(&stop);
    
    // Benchmark coalesced
    cudaEventRecord(start);
    for (int i = 0; i < 100; i++) {
        coalesced<<<gridSize, blockSize>>>(d_output, d_input, N);
    }
    cudaEventRecord(stop);
    cudaEventSynchronize(stop);
    
    float ms1 = 0;
    cudaEventElapsedTime(&ms1, start, stop);
    
    // Benchmark non-coalesced
    cudaEventRecord(start);
    for (int i = 0; i < 100; i++) {
        nonCoalesced<<<gridSize, blockSize>>>(d_output, d_input, N);
    }
    cudaEventRecord(stop);
    cudaEventSynchronize(stop);
    
    float ms2 = 0;
    cudaEventElapsedTime(&ms2, start, stop);
    
    printf("=== Memory Access Pattern Benchmark ===\n");
    printf("Coalesced:     %.3f ms\n", ms1);
    printf("Non-coalesced: %.3f ms\n", ms2);
    printf("Slowdown:      %.2fx\n", ms2 / ms1);
    
    cudaFree(d_input);
    cudaFree(d_output);
    free(h_input);
    free(h_output);
    cudaEventDestroy(start);
    cudaEventDestroy(stop);
}

int main() {
    benchmarkMemoryAccess(10000000);
    return 0;
}

4.2 共享内存矩阵乘法

cpp

#include <cuda_runtime.h>
#include <stdio.h>

#define TILE_SIZE 16

// Naive matrix multiplication (global memory only)
__global__ void matMulNaive(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;
    }
}

// Optimized matrix multiplication using shared memory
__global__ void matMulOptimized(const float* A, const float* B, float* C,
                                 int M, int N, int K) {
    __shared__ float As[TILE_SIZE][TILE_SIZE];
    __shared__ float Bs[TILE_SIZE][TILE_SIZE];
    
    int tx = threadIdx.x;
    int ty = threadIdx.y;
    int row = blockIdx.y * TILE_SIZE + ty;
    int col = blockIdx.x * TILE_SIZE + tx;
    
    float sum = 0.0f;
    
    int numTiles = (K + TILE_SIZE - 1) / TILE_SIZE;
    
    for (int t = 0; t < numTiles; t++) {
        // Load A tile
        int aCol = t * TILE_SIZE + tx;
        if (row < M && aCol < K) {
            As[ty][tx] = A[row * K + aCol];
        } else {
            As[ty][tx] = 0.0f;
        }
        
        // Load B tile
        int bRow = t * TILE_SIZE + ty;
        if (bRow < K && col < N) {
            Bs[ty][tx] = B[bRow * N + col];
        } else {
            Bs[ty][tx] = 0.0f;
        }
        
        __syncthreads();
        
        // Compute partial sum
        for (int k = 0; k < TILE_SIZE; k++) {
            sum += As[ty][k] * Bs[k][tx];
        }
        
        __syncthreads();
    }
    
    if (row < M && col < N) {
        C[row * N + col] = sum;
    }
}

// Host function to benchmark
void benchmarkMatMul(int M, int N, int K) {
    size_t sizeA = M * K * sizeof(float);
    size_t sizeB = K * N * sizeof(float);
    size_t sizeC = M * N * sizeof(float);
    
    float *h_A = (float*)malloc(sizeA);
    float *h_B = (float*)malloc(sizeB);
    float *h_C = (float*)malloc(sizeC);
    
    for (int i = 0; i < M * K; i++) h_A[i] = (float)i / (M * K);
    for (int i = 0; i < K * N; i++) h_B[i] = (float)i / (K * N);
    
    float *d_A, *d_B, *d_C;
    cudaMalloc((void**)&d_A, sizeA);
    cudaMalloc((void**)&d_B, sizeB);
    cudaMalloc((void**)&d_C, sizeC);
    
    cudaMemcpy(d_A, h_A, sizeA, cudaMemcpyHostToDevice);
    cudaMemcpy(d_B, h_B, sizeB, cudaMemcpyHostToDevice);
    
    dim3 blockSize(TILE_SIZE, TILE_SIZE);
    dim3 gridSize((N + TILE_SIZE - 1) / TILE_SIZE, 
                  (M + TILE_SIZE - 1) / TILE_SIZE);
    
    cudaEvent_t start, stop;
    cudaEventCreate(&start);
    cudaEventCreate(&stop);
    
    // Warm-up
    matMulOptimized<<<gridSize, blockSize>>>(d_A, d_B, d_C, M, N, K);
    cudaDeviceSynchronize();
    
    // Benchmark optimized
    cudaEventRecord(start);
    for (int i = 0; i < 50; i++) {
        matMulOptimized<<<gridSize, blockSize>>>(d_A, d_B, d_C, M, N, K);
    }
    cudaEventRecord(stop);
    cudaEventSynchronize(stop);
    
    float ms1 = 0;
    cudaEventElapsedTime(&ms1, start, stop);
    
    // Benchmark naive
    cudaEventRecord(start);
    for (int i = 0; i < 50; i++) {
        matMulNaive<<<gridSize, blockSize>>>(d_A, d_B, d_C, M, N, K);
    }
    cudaEventRecord(stop);
    cudaEventSynchronize(stop);
    
    float ms2 = 0;
    cudaEventElapsedTime(&ms2, start, stop);
    
    printf("=== Matrix Multiplication Benchmark (%dx%dx%d) ===\n", M, N, K);
    printf("Optimized (shared mem): %.3f ms\n", ms1 / 50);
    printf("Naive (global mem):     %.3f ms\n", ms2 / 50);
    printf("Speedup:                %.2fx\n", ms2 / ms1);
    
    cudaMemcpy(h_C, d_C, sizeC, cudaMemcpyDeviceToHost);
    
    cudaFree(d_A);
    cudaFree(d_B);
    cudaFree(d_C);
    free(h_A);
    free(h_B);
    free(h_C);
    cudaEventDestroy(start);
    cudaEventDestroy(stop);
}

int main() {
    benchmarkMatMul(1024, 1024, 1024);
    return 0;
}

4.3 常量内存与纹理内存示例

cpp

#include <cuda_runtime.h>
#include <stdio.h>

// Constant memory declaration
__constant__ float gaussianFilter[9];

// Kernel using constant memory
__global__ void constantMemKernel(const float* input, float* output, 
                                   int width, int height) {
    int x = blockIdx.x * blockDim.x + threadIdx.x;
    int y = blockIdx.y * blockDim.y + threadIdx.y;
    
    if (x < width && y < height) {
        float sum = 0.0f;
        int filterRadius = 1;
        
        for (int fy = -filterRadius; fy <= filterRadius; fy++) {
            for (int fx = -filterRadius; fx <= filterRadius; fx++) {
                int px = min(width - 1, max(0, x + fx));
                int py = min(height - 1, max(0, y + fy));
                
                int filterIdx = (fy + filterRadius) * 3 + (fx + filterRadius);
                sum += input[py * width + px] * gaussianFilter[filterIdx];
            }
        }
        
        output[y * width + x] = sum;
    }
}

// Texture reference for texture memory
texture<float, 2, cudaReadModeElementType> texRef;

// Kernel using texture memory
__global__ void textureMemKernel(float* output, int width, int height) {
    int x = blockIdx.x * blockDim.x + threadIdx.x;
    int y = blockIdx.y * blockDim.y + threadIdx.y;
    
    if (x < width && y < height) {
        // Texture fetch with automatic caching and spatial locality
        float val = tex2D(texRef, (float)x, (float)y);
        output[y * width + x] = val * 2.0f;
    }
}

// Host setup for constant memory
void setupConstantMemory() {
    float h_filter[9] = {
        1.0f/16, 2.0f/16, 1.0f/16,
        2.0f/16, 4.0f/16, 2.0f/16,
        1.0f/16, 2.0f/16, 1.0f/16
    };
    cudaMemcpyToSymbol(gaussianFilter, h_filter, 9 * sizeof(float));
}

// Host setup for texture memory
void setupTextureMemory(const float* d_data, int width, int height) {
    cudaChannelFormatDesc channelDesc = cudaCreateChannelDesc<float>();
    cudaArray* cuArray;
    cudaMallocArray(&cuArray, &channelDesc, width, height);
    
    cudaMemcpyToArray(cuArray, 0, 0, d_data, 
                      width * height * sizeof(float), 
                      cudaMemcpyDeviceToDevice);
    
    texRef.addressMode[0] = cudaAddressModeClamp;
    texRef.addressMode[1] = cudaAddressModeClamp;
    texRef.filterMode = cudaFilterModePoint;
    texRef.normalizedCoords = false;
    
    cudaBindTextureToArray(texRef, cuArray, channelDesc);
}

int main() {
    int width = 512;
    int height = 512;
    size_t size = width * height * sizeof(float);
    
    float *h_data = (float*)malloc(size);
    for (int i = 0; i < width * height; i++) {
        h_data[i] = (float)(i % 256);
    }
    
    float *d_data, *d_output;
    cudaMalloc((void**)&d_data, size);
    cudaMalloc((void**)&d_output, size);
    
    cudaMemcpy(d_data, h_data, size, cudaMemcpyHostToDevice);
    
    // Setup constant memory
    setupConstantMemory();
    
    // Run constant memory kernel
    dim3 blockSize(16, 16);
    dim3 gridSize((width + 15) / 16, (height + 15) / 16);
    
    constantMemKernel<<<gridSize, blockSize>>>(d_data, d_output, width, height);
    cudaDeviceSynchronize();
    
    printf("Constant memory kernel completed.\n");
    
    // Setup and run texture memory kernel
    setupTextureMemory(d_data, width, height);
    
    textureMemKernel<<<gridSize, blockSize>>>(d_output, width, height);
    cudaDeviceSynchronize();
    
    printf("Texture memory kernel completed.\n");
    
    cudaFree(d_data);
    cudaFree(d_output);
    free(h_data);
    
    return 0;
}

五、深入思考

5.1 历史背景

GPU内存架构经历了显著的演进。早期的GPU只有全局显存和寄存器，共享内存在Fermi架构中才真正成熟。常量内存和纹理内存在Tesla架构中引入，最初为图形渲染设计，后来发现它们在通用计算中也有广泛应用。

L1缓存在Fermi架构中首次引入，之前GPU完全依赖程序员管理的共享内存来实现数据复用。L2缓存在Fermi中为768KB，在后续架构中不断扩大，Ampere A100达到40MB。

统一内存在CUDA 6.0中引入，最初基于简单的页面映射机制。Pascal架构引入了页面迁移引擎，Volta进一步改进了缺页异常处理，使得统一内存在生产环境中更加实用。

5.2 前沿方向

HBM3与更高带宽：新一代HBM3显存提供超过3 TB/s的带宽，将进一步缓解内存瓶颈。但这也要求程序员更加注重合并访问，以充分利用高带宽。

片上内存扩展：未来的GPU可能会增加共享内存和L1缓存的容量，减少对外部显存的依赖。3D堆叠技术使得在GPU芯片上集成更多SRAM成为可能。

智能预取：硬件预取器正在变得更加智能，可以识别复杂的访问模式并提前加载数据。程序员可能需要考虑如何编写预取友好的代码。

持久化线程与内存：CUDA合作组（Cooperative Groups）引入了持久化线程的概念，线程可以长时间驻留，维护状态。这对内存管理提出了新的挑战和机遇。

5.3 延伸阅读

NVIDIA CUDA C Programming Guide - Memory Model

"Inside the NVIDIA Volta GPU Architecture" - IEEE Micro

NVIDIA Developer Blog: "How to Access Global Memory Efficiently in CUDA"

CUDA Best Practices Guide - Memory Optimizations

六、实验笔记

实验1：合并访问性能量化

目的：测量不同访问模式的带宽利用率。

步骤：

实现顺序访问、步长访问、随机访问三种kernel

使用Nsight Compute测量每种模式的DRAM利用率

计算理论带宽与实际带宽的比值

预期：顺序访问应达到80%以上利用率，步长访问随步长增大而降低，随机访问最低。

实验2：共享内存vs全局内存性能对比

目的：量化共享内存优化的收益。

步骤：

实现两个版本的矩阵乘法：一个仅使用全局内存，一个使用共享内存优化

测量不同矩阵大小下的执行时间

分析共享内存减少的全局内存访问次数

预期：对于大矩阵，共享内存版本应有2-5倍加速。小矩阵时差异较小，因为数据可能已在缓存中。

七、本章小结

本章全面介绍了GPU的内存层次结构，从最快的寄存器到最慢的全局显存，每种内存类型都有其独特的特性和适用场景。我们深入分析了全局内存的合并访问机制，这是GPU内存优化的核心：相邻线程访问相邻地址可以实现最高的带宽利用率。

共享内存是GPU编程中最强大的优化工具，它允许程序员显式管理数据在SM内的存储和复用。通过共享内存，可以实现高效的线程间通信和数据复用，大幅减少全局内存访问。但共享内存的使用需要注意bank conflict问题，合理的填充和数据重排可以避免这一陷阱。

常量内存和纹理内存是针对特定访问模式优化的专用存储器。常量内存的广播机制使得warp内的一致读取极其高效；纹理内存的空间局部性优化使其在图像处理和卷积应用中表现出色。

L1/L2缓存由硬件管理，为程序员提供了透明的加速。虽然无法直接控制缓存内容，但通过编写缓存友好的代码（良好的空间和时间局部性），可以间接提高缓存命中率。

统一内存简化了编程模型，但页面迁移开销可能影响性能。对于性能敏感的应用，显式的内存管理仍然是更好的选择。

通过三个代码示例，我们实践了合并访问对比、共享内存矩阵乘法和常量/纹理内存的使用，将理论知识转化为实际的编程技能。

八、习题

合并访问分析：对于float数组，thread i访问data[i * 3]。这是否是合并访问？如果不是，带宽利用率大约是多少？如何重新组织数据以实现合并访问？

共享内存容量计算：假设SM有192KB的片上内存，配置为112KB共享内存 + 80KB L1缓存。如果每个Block使用4KB共享内存，SM最多可以驻留多少个Block？如果改为64KB共享内存 + 128KB L1缓存，结果如何？

Bank Conflict分析：共享内存数组float shared[32][32]，thread (tx, ty)访问shared[tx][ty]。是否存在bank conflict？如果改为shared[ty][tx]呢？如果数组改为float shared[32][33]，再次访问shared[ty][tx]，情况如何？

常量内存适用性：以下场景是否适合使用常量内存？说明理由。

a) 存储1000个元素的查找表，所有线程访问相同索引

b) 存储1000个元素的查找表，每个线程访问不同索引

c) 存储kernel的配置参数（10个float）

缓存命中率估算：一个kernel遍历一个10MB的数组，每个元素访问一次。L2缓存为8MB。估算L2缓存命中率。如果数组大小为5MB呢？

矩阵转置优化：编写一个使用共享内存优化的矩阵转置kernel，避免bank conflict。与非优化版本比较性能。

统一内存页面迁移：编写一个测试程序，交替在主机和设备上访问统一内存分配的数组。测量页面迁移的次数和开销。使用cudaPrefetchAsync优化后，性能改善多少？

纹理内存应用：除了图像处理，纹理内存在哪些通用计算场景中可能有用？考虑具有空间局部性的数据访问模式。

内存带宽计算：Ampere A100的显存带宽为1.6 TB/s。如果一个kernel每次内存事务传输128字节，每秒最多可以进行多少次内存事务？如果warp内有32个线程，每个线程请求4字节，理想情况下需要多少次内存事务？

综合优化：设计一个高效的并行直方图计算kernel。说明如何使用共享内存减少全局内存原子操作，如何处理bin冲突，以及如何优化内存访问模式。