目录
一：概述
二：内存访问效率的重要性
三：RoofLine 模型
四：CUDA 内存类型
五：用分块技术减少全局内存访问
六：分块矩阵乘法
七：边界处理
八：内存使用对占用率的影响（occupancy）
九：总结
一：概述
到目前为止，我们已经学会了如何编写 CUDA 核函数，以及如何设置和分配大量线程来执行核函数。我们还了解了当前 GPU 硬件的计算架构，以及线程在硬件上调度执行过程。在本章中，我们将重点关注 GPU 的片上(on-chip)内存架构，并研究如何组织和存放数据，以便这些数据能够被大量线程高效的访问。
到目前为止，我们所学习的 CUDA 核函数可能只达到底层硬件性能的一小部分。之所以性能不佳，是因为通常使用片外（off-chip）内存，即全局内存。全局内存往往具有较长的访问延迟（数百个时钟周期）和有限的访问带宽。尽管理论上当有许多线程执行时可以容许有较长的内存访问延迟。 但很容易出现这样的情况：当大量线程同时去访问全局内存时，由于全局内存带宽有限，容易产生交通阻塞，导致除极少数线程外，其他所有线程都无法及时访问内存。 这样，流式多处理器（SM）中的部分计算核（cores）就会空闲。为了避免这种情况，GPU 提供了大量的片上（on-chip）内存资源，用于数据的访问，从而消除了访问全局内存的交通阻塞情况。在本章中，我们将研究如何使用不同的内存类型来提高 CUDA 核函数的执行性能。
二：内存访问效率的重要性
回顾一下前一章矩阵乘法的例子（CUDA编程04 - GPU计算架构和线程调度）， 我们可以在矩阵乘法的核函数中找出矩阵乘法执行最多的那部分代码，通过计算这部分代码的预期性能来说明内存访问对效率的影响。 下图中复制了这部分代码。就执行时间而言，核函数中最重要的部分是在 for 循环中执行 M 行与 N 列点积的代码。
在循环的每次迭代中，会进行两次全局内存访问，以执行一次浮点乘法和一次浮点加法。全局内存访问从 M 和 N 数组中获取元素。浮点乘法操作将这两个元素相乘，而浮点加法操作将乘积累加到 Pvalue 中。因此，浮点操作（FLOP）与从全局内存访问的字节（B）之比为 2 FLOP / 8 B，或 0.25 FLOP/B。我们将这个比率称为计算与全局内存访问比率，定义为在程序的某个区域内每次从全局内存访问的字节所执行的 FLOP 数量。在文献中，这个比率有时也被称为算术强度或计算强度。
计算与全局内存访问比率对CUDA核函数的性能有重大影响。例如，Ampere A100 GPU的峰值全局内存带宽为1555 GB/秒。由于矩阵乘法内核执行0.25 OP/B，全局内存带宽将核函数可以执行的单精度FLOP吞吐量限制为每秒389 giga FLOPS（GFLOPS），这是通过将1555 GB/秒乘以0.25 FLOP/B获得的。然而，389 GFLOPS仅占A100 GPU峰值单精度操作吞吐量的2%，即19,500 GFLOPS。A100还配备了称为张量核心的专用单元，适用于加速矩阵乘法操作。如果考虑A100的张量核心峰值单精度浮点吞吐量为156,000 GFLOPS，389 GFLOPS仅占峰值的0.25%。因此，矩阵乘法内核的执行受到从内存到GPU核心的数据传输速率的严重限制。我们将执行速度受内存带宽限制的程序称为内存瓶颈程序。
为了提高该核函数的性能，我们需要通过减少全局内存访问次数来提高内核的计算与全局内存访问比率。例如，要充分利用 A100 GPU 提供的 19,500GFLOPS 性能，至少需要（19,500 GOP/秒）/（1555 GB/秒）=12.5 OP/B。这一比率意味着，每访问一个 4 字节浮点数值，必须执行约 50 次浮点运算！能在多大程度上实现这一比率取决于当前计算中的内在数据重用情况。读者可参阅 “Roofline 模型 ”，了解分析程序潜在性能与其计算强度关系的有用模型。
我们将看到，矩阵乘法提供了减少全局内存访问的机会，而这些机会可以通过相对简单的技术来把握。矩阵乘法函数的执行速度可能会有数量级的变化，这取决于减少全局内存访问的程度。本章将介绍一种减少全局内存访问次数的常用技术，并在矩阵乘法中演示该技术。
三：RoofLine 模型
RoofLine模型是一种可视化模型，用于评估应用程序相对于其运行的硬件极限的性能。RoofLine模型的基本示例如下。
在 x 轴上，我们用 FLOP/B 来衡量算术或计算强度。它反映了应用程序每加载一个字节数据所做的工作量。在 Y 轴上，计算性能以 GFLOPS 为单位。图中的两条线反映了硬件极限。水平线表示硬件计算性能（GFLOPS）的极限。从原点开始的斜线表示硬件内存带宽性能的极限。图中的“点”代表一个应用，X 轴为其运行强度，Y 轴为其达到的计算性能。当然，这些点会在两条线之下，因为它们无法达到比硬件峰值更高的性能。
从点与两条线的相对位置可以看出应用程序的效率。接近两条线的点表示应用程序正在高效使用内存带宽或计算单元，而远低于两条线的应用程序则表示资源使用效率低下。这两条线的交点代表计算强度值，在此值上，应用程序从内存带宽瓶颈过渡到计算瓶颈。计算强度较低的应用程序受内存带宽限制，无法达到峰值性能。计算强度较高的应用属于计算瓶颈，不受内存带宽限制。
例如，A1 和 A2 点都代表内存带宽瓶颈的应用程序，而 A3 则代表计算瓶颈的应用程序。A1 能有效利用资源，并在接近内存带宽峰值的情况下运行，而 A2 则不然。对 A2 来说，可能还有额外优化的余地，可以通过提高内存带宽利用率来提高性能，但对 A1 来说，提高性能的唯一方法就是提高应用程序的计算强度。
四：CUDA 内存类型
CUDA设备包含几种类型的内存，可以帮助程序员提高计算与全局内存访问的比率。下图显示了这些CUDA设备内存。在图的底部，我们可以看到全局内存和常量内存。这两种类型的内存都可以被主机写入（W）和读取（R）。全局内存也可以被设备写入和读取，而常量内存支持设备的短延迟、高带宽只读访问。我们在后面文章中将详细讨论常量内存。
另一种类型的内存是局部内存，它也可以被读取和写入。局部内存实际上放置在全局内存中，具有类似的访问延迟，但它不在线程之间共享。每个线程都有自己的一部分全局内存，作为其私有局部内存，用于存放线程私有但无法在寄存器中分配的数据。这些数据包括静态分配的数组、寄存器装不下的变量以及线程调用栈上的变量。
上图中的寄存器和共享内存是片上存储器。这些类型的内存中驻留的变量可以以非常高的速度以高度并行的方式访问。寄存器分配给单独的线程；每个线程只能访问自己的寄存器（请参阅“CPU与GPU寄存器架构”）。内核函数通常使用寄存器来保存每个线程私有的频繁访问的变量。共享内存分配给线程块；块中的所有线程都可以访问为该块声明的共享内存变量。共享内存是线程通过共享输入数据和中间结果进行协作的有效手段。通过在CUDA内存类型之一中声明CUDA变量，CUDA程序员决定了变量的可见性和访问速度。
CPU和GPU之间不同的设计目标导致了不同的寄存器架构。正如我们在前一篇文章中看到的，当CPU在不同线程之间进行上下文切换时，它们将即将退出线程的寄存器保存到内存中，并从内存中恢复即将进入线程的寄存器。相比之下，GPU通过将所有调度在处理块上的线程的寄存器保存在处理块的寄存器文件中，实现了零开销调度。因此，线程组之间的切换是瞬时的，因为即将进入线程的寄存器已经在寄存器文件中。因此，GPU寄存器文件需要比CPU寄存器文件大得多。
我们在前一篇文章中还看到，GPU支持动态资源分配，其中一个SM可以为每个线程分配少量寄存器并执行大量线程，或者为每个线程分配更多寄存器并执行较少线程。因此，GPU寄存器文件需要设计成支持这种寄存器的动态分配。相比之下，CPU寄存器架构为每个线程分配固定数量的寄存器，而不考虑线程对寄存器的实际需求。
为了充分理解寄存器、共享内存和全局内存之间的差异，我们需要更详细地讨论这些不同内存类型在现代处理器中的实现和使用。如我们在前一篇文章“Warps与SIMD硬件”中讨论的，几乎所有现代处理器都源于约翰·冯·诺依曼在1945年提出的模型，如下图所示。CUDA设备也不例外。CUDA设备中的全局内存映射到下图中的内存（Memory）框。处理器（Processing Unit）框对应于我们今天通常看到的处理器芯片边界。全局内存位于处理器芯片之外，并采用DRAM技术实现，这意味着访问延迟较长，访问带宽相对较低。寄存器对应于冯·诺依曼模型的“寄存器文件”。寄存器文件位于处理器芯片上，这意味着与全局内存相比，访问延迟非常短，访问带宽显著更高。在典型设备中，所有SM的寄存器文件的聚合访问带宽至少比全局内存高两个数量级。此外，每当一个变量存储在寄存器中时，其访问将不再消耗off-chip的全局内存带宽。这将反映为计算与访存比值增加。
更微妙的一点是，每次访问寄存器的指令比对全局内存访问的指令少。大多数现代处理器中的算术指令都有“内置”的寄存器操作数。例如，浮点加法可以采用以下形式:
其中 r2 和 r3 是寄存器编号，指定寄存器文件中输入操作数的位置。。存储浮点加法结果值的存放位置由 r1 指定。因此，当算术指令的操作数位于寄存器中时，不需要额外的指令就能将操作数值提供给算术和逻辑单元 (ALU)，算术计算就已经可以做完了。
同时，如果一个操作数值在全局内存中，处理器需要把数据先加载到寄存器中，使操作数对ALU可用。例如，如果浮点加法指令的第一个操作数在全局内存，所涉及的指令可能看起来像下面的例子:
加载指令将内存偏移量放在 r4 的寄存将器中，以形成操作数值的地址。然后，它访问全局内存并将值放入寄存器 r2。一旦操作数值在 r2 中，fadd 指令使用 r2 和 r3 中的值执行浮点加法，并将结果放入 r1。由于处理器每个时钟周期只能获取和执行有限数量的指令，带有额外加载到寄存器的处理时间可能会比没有的更长。这是将操作数放入寄存器可以提高执行速度的另一个原因。
最后，将操作数放在寄存器中还有一个微妙的原因。 在现代计算机中，从寄存器文件访问一个值所消耗的能量至少比从全局存储器访问一个值所消耗的能量至少低一个数量级。 从寄存器中访问数值比从全局存储器中访问数值在能效方面具有极大的优势。 我们将在现代计算机看到，访问这两种硬件结构的速度和能耗差异的更多细节。 另一方面，正如我们即将看到的，在当今的 GPU 中，每个线程可使用的寄存器数量非常有限。 正如我们在前一篇文章中看到的，如果能减少线程的占用率(occupancy)，就能提高应用程序的性能。如果full-occupancy的情况下寄存器使用量超过了硬件最大寄存器数量，应用程序的占用率就会降低。 因此，我们还需要尽可能避免超量使用这一有限资源。
下图显示了 CUDA 设备中的共享内存和寄存器。 虽然都是片上存储器，但它们在功能和访问成本上有很大不同。 共享内存被设计为驻留在处理器芯片上的内存空间的一部分。当处理器访问共享内存中的数据时，需要执行内存加载操作，就像访问全局内存中的数据一样。不过，由于共享内存位于芯片上，因此它的访问延迟和吞吐量要比全局存储器低得多。 由于需要执行加载操作，共享内存的延迟更长，带宽更低。 在计算机体系结构术语中，共享内存是一种折衷内存类型。
在CUDA中，共享内存和寄存器之间一个重要的区别是，存储在共享内存中的变量可以被块中的所有线程访问。这与寄存器数据相对，后者是线程私有的。也就是说，共享内存旨在支持块内线程之间高效、高带宽的数据共享。如上图所示，CUDA设备的SM通常使用多个处理单元，以允许多个线程同时进行访问。块中的线程可以分布在这些处理单元上。因此，这些CUDA设备中的共享内存的硬件实现通常被设计为允许多个处理单元同时访问其内容，以支持块内线程之间的高效数据共享。我们将学习几种重要类型的并行算法，这些算法可以从线程之间的高效数据共享中获得巨大收益。
现在应该清楚的是，寄存器、本地内存、共享内存和全局内存具有不同的功能、延迟和带宽。因此，了解如何声明变量以使预期内存类型是很重要的。下图展示了将程序变量声明为各种内存类型的CUDA语法。每个这样的声明还为其声明的CUDA变量赋予了作用域和生命周期。作用域确定可以访问该变量的线程集：仅单个线程、一个块的所有线程或所有网格的所有线程。如果变量的作用域是单个线程，则每个线程将为该变量创建一个私有版本；每个线程只能访问其私有版本的变量。例如，如果一个核函数声明了一个作用域为线程的变量，并且它以一百万个线程启动，那么将创建一百万个版本的变量，以便每个线程初始化和使用其自己的版本的变量。
生命周期指的是程序执行期间变量可用的部分：要么在网格的生命周期内，要么在整个应用程序生命周期内。如果变量的生命周期在网格的生命周期内，则必须在核函数体内声明，并且仅可由内核的代码使用。如果内核被多次调用，则变量的值不会在这些调用之间持续保存。每次调用必须初始化变量以便使用。另一方面，如果变量的生命周期贯穿整个应用程序，则必须在任何函数体外声明。这些变量的内容在应用程序执行期间持续保存，并且对所有核函数可用。
我们将非数组的变量称为标量变量。如上图所示，所有在核函数和设备函数中声明的自动标量变量都被放置在寄存器中。这些自动变量的作用域在各个线程内。当核函数声明一个自动变量时，为执行该核函数的每个线程生成一个该变量的私有副本。当线程终止时，所有其自动变量都不复存在。在之前图中的变量blurRow、blurCol、curRow、curCol、pixels和pixVal都是自动变量，属于这一类别。请注意，访问这些变量的速度非常快并且是并行的，但必须小心不要超过硬件实现中寄存器存储的有限容量。使用大量寄存器可能会对每个SM的占用率产生负面影响，正如我们在前面一篇文章CUDA编程04 - GPU计算架构和线程调度中所看到的。
标量数组变量不存储在寄存器中 ，而是存储在线程的本地内存中，可能会导致较长的访问延迟和潜在的访问拥塞。 与标量变量一样，这些数组的范围也仅限于单个线程。 也就是说，每个标量数组都是为每个线程创建并使用的。 一旦线程终止执行，其标量数组变量的内容就不复存在。 根据我们的经验，核函数和设备函数中很少需要使用标量数组变量。
如果变量声明前面有 __shared__ 关键字（每个“__”由两个“_”字符组成），则在 CUDA 中声明了一个共享变量。在声明中也可以在 __shared__ 前添加可选的 __device__ 以达到相同的效果。这种声明通常在核函数或设备函数内进行。共享变量驻留在共享内存中。共享变量的作用域在一个线程块内；也就是说，块中的所有线程都看到同一个版本的共享变量。在核函数执行期间，为每个块创建并使用共享变量的私有版本。共享变量的生命周期是在核函数执行期间。当一个核函数终止其网格的执行时，其共享变量的内容将不复存在。正如我们之前讨论的，共享变量是块内线程相互协作的有效手段。从共享内存访问共享变量是极其快速和高度并行的。CUDA 程序员通常使用共享变量来保存在内核执行阶段中频繁使用和重用的全局内存数据的一部分。正如我们将在后面通过矩阵乘法演示的那样，可能需要调整用于创建执行分阶段的算法，利用共享内存，以每次计算全局内存数据的一小部分。
如果变量声明前面有关键字 __constant__（每个“__”由两个“_”字符组成），则在CUDA中声明一个常量变量。也可以在 __constant__ 前添加一个可选的 __device__ 以实现相同的效果。常量变量的声明必须在任何函数体外。常量变量的作用域是所有网格，这意味着所有网格中的所有线程都看到相同版本的常量变量。常量变量的生命周期是整个应用程序的执行。常量变量通常用于提供输入值给核函数的变量。常量变量的值不能被内核函数代码更改。常量变量存储在全局内存中，但会被缓存以实现高效访问。通过适当的访问模式，访问常量内存是非常快速和并行的。目前，应用程序中常量变量的总大小限制为65,536字节。可能需要将输入数据量拆分以适应此限制。我们将在后面的卷积计算中演示常量内存的使用。
一个仅由关键字 __device__（每个“__”由两个“_”字符组成）前置声明的变量是全局变量，将被放置在全局内存中。访问全局变量的速度较慢。随着更近期设备中缓存的引入，访问全局变量的延迟和吞吐量得到了改善。全局变量的一个重要优点是它们对所有核函数的所有线程都是可见的。它们的内容在整个执行过程中也会持续存在。因此，全局变量可以作为线程在不同块之间协作的手段。然而，必须注意的是，目前没有简单的方法可以在来自不同线程块的线程之间进行同步，或确保线程在访问全局内存时的数据一致性，除了使用原子操作或终止当前内核执行。因此，全局变量通常用于在一个核函数调用与另一个核函数调用之间传递信息。
在CUDA中，指针可以用于指向全局内存中的数据对象。指针使用在内核和设备函数中出现的两种典型方式。首先，如果对象由主机函数分配，则对象的指针通过内存分配API函数（如cudaMalloc）初始化，并可以作为参数传递给核函数，正如我们在前面CUDA编程02 - 数据并行介绍和CUDA编程03 - 多维数据并行中看到的那样。第二种使用方式是将声明在全局内存中的变量的地址赋值给指针变量。例如，在内核函数中的语句{float *ptr = &GlobalVar;}将GlobalVar的地址赋给自动指针变量ptr。读者应参考CUDA编程指南，以了解在其他内存类型中使用指针的方法。
五：用分块技术减少全局内存访问
我们在使用CUDA设备内存时面临着内在的权衡：全局内存容量大但速度慢，而共享内存容量小但速度快。一种常见的策略是将数据分割成称为“瓷砖”的子集，以便每个瓷砖（tile）都能放入共享内存中。“瓷砖”这个术语借鉴了这样的类比：一个大的墙面（即全局内存数据）可以用小的瓷砖（即可以放入共享内存中的子集）覆盖。一个重要的标准是，这些瓷砖上的核函数计算可以相互独立地进行。请注意，对于任意的核函数，并非所有数据结构都可以分割成瓷砖（tile），以下分割“瓷砖”简称分块。
分块的概念可以通过 CUDA编程03 - 多维数据并行 中的矩阵乘法示例来说明。为了方便参考，我们在下图中复制了该示例。为简洁起见，我们将P[y*Width +x]、M[y*Width +x]和N[y*Width +x]简写为P[y,x]、M[y,x]和N[y,x]。该示例假设我们使用4x2的块来计算P矩阵。P矩阵中的黑色框定义了由每个块处理的P元素。下图突出显示了块[0,0]的四个线程所进行的计算。这四个线程计算P[0,0]、P[0,1]、P[1,0]和P[1,1]。线程[0,0]和线程[0,1]对块[0,0]的M和N元素的访问用黑色箭头突出显示。例如，线程[0,0]读取M[0,0]和N[0,0]，接着是M[0,1]和N[1,0]，然后是M[0,2]和N[2,0]，最后是M[0,3]和N[3,0]。
图5.6显示了由块[0，0]中的所有线程对全局内存的访问。线程沿着垂直方向从上到下，访问时间沿水平方向上从左到右增加。注意，每个线程在其执行期间访问M的四个元素和N的四个元素。在突出显示的四个线程中，在它们访问的m和n元素中有显著的重叠。例如，线程[0，1]和线程[1，1]都访问N[0，1]以及N的第1列的其余部分。
CUDA编程03 - 多维数据并行 矩阵乘法核函数的编写方式使得thread[0,0]和thread[0,1]都从全局内存访问M的第0行元素。如果我们能够以某种方式让thread[0,0]和thread[0,1]协作，使得这些M元素仅从全局内存加载一次，我们就可以将对全局内存的总访问次数减少一半。事实上，我们可以看到在执行block[0,0]期间，每个M和N元素正好被访问两次。因此，如果我们能够让所有四个线程在访问全局内存时进行协作，我们可以将对全局内存的流量减少一半。
读者应该验证得到在矩阵乘法示例中，全局内存访问的减少与所使用的块的维度成正比。使用宽度为Width*Width的块时，全局内存访问的减少为Width。也就是说，如果我们使用16 x 16的块，通过线程之间的协作，我们可以将全局内存流量减少到原始水平的1/16。
我们现在介绍一种分块矩阵乘法算法。基本思路是让线程协作地将 M 和 N 元素的子集加载到共享内存中，然后再单独使用这些元素进行点积计算。请记住，共享内存的大小非常小，因此在将这些 M 和 N 元素加载到共享内存时，必须小心不要超过共享内存的容量。这可以通过将 M 和 N 矩阵划分为更小的块来实现。这些块的大小选择得足以适应共享内存。在最简单的形式中，块的维度等于分块的维度，如下图所示。
在上图中，我们将M和N划分为2 x 2的瓦片，如粗线所示。每个线程执行的点积计算现在被分为多个阶段。在每个阶段，块中的所有线程协作将M的一个瓦片和N的一个瓦片加载到共享内存中。这可以通过让块中的每个线程将一个M元素和一个N元素加载到共享内存中来完成，如下图所示。下图的每一行显示了一个线程的执行活动。请注意，时间是从左到右推进的。我们只需要显示块0,0中线程的活动；其他块的行为都是相同的。M元素的共享内存数组称为Mds。N元素的共享内存数组称为Nds。在阶段1开始时，块0,0中的四个线程协作将一个M瓦片加载到共享内存中：线程0,0将M0,0加载到Mds0,0，线程0,1将M0,1加载到Mds0,1，线程1,0将M1,0加载到Mds1,0，线程1,1将M1,1加载到Mds1,1。这些加载在图5.8的第二列中显示。N的一个瓦片也以类似的方式加载，如下图的第三列所示。
在将 M 和 N 的两个块加载到共享内存后，这些元素用于计算点积。请注意，共享内存中的每个值被使用两次。例如，由线程1,1 加载到 Mds1,1 的 M1,1 值被使用两次，一次由线程1,0，另一次由线程1,1。通过将每个全局内存值加载到共享内存中，以便可以多次使用，我们减少了对全局内存的访问次数。在这种情况下，我们将对全局内存的访问次数减少了一半。读者应验证，如果块是 N x N 元素，则减少的次数是 N 倍。
请注意，每个点积的计算现在分为两个阶段进行，如图5.8中的阶段1和阶段2所示。在每个阶段中，每个线程将输入矩阵元素的两个对的乘积累加到Pvalue变量中。请注意，Pvalue是一个自动变量，因此为每个线程生成一个私有版本。我们添加了下标只是为了澄清这些是为每个线程创建的不同实例的Pvalue变量。第一阶段的计算在图5.8的第四列中显示，第二阶段在第七列中显示。一般来说，如果输入矩阵的维度为Width，瓷砖大小为TILE_WIDTH，则点积将在Width/TILE_WIDTH个阶段中进行。这些阶段的创建是减少对全局内存访问的关键。每个阶段专注于输入矩阵值的一个小子集，线程可以协作将该子集加载到共享内存中，并使用共享内存中的值来满足它们在该阶段的重叠输入需求。
请注意，Mds 和 Nds 在各个阶段之间是重用的。在每个阶段，使用相同的 Mds 和 Nds 来保存该阶段使用的 M 和 N 元素的子集。这允许一个更小的共享内存来处理大多数对全局内存的访问。这是因为每个阶段专注于输入矩阵元素的一个小子集。这种集中访问行为被称为局部性。当算法表现出局部性时，就有机会使用小型高速内存来处理大多数访问，并将这些访问从全局内存中移除。局部性在实现多核 CPU 和多线程 GPU 的高性能方面同样重要。我们将在下一篇文章中回到局部性的概念。
六：分块矩阵乘法
现在我们介绍一个使用共享内存做矩阵乘法的核函数的例子，来减少对全局内存的访问。核函数下图所示。在图中，第04行和第05行声明了Mds和Nds，分别为共享内存数组。共享内存的作用域是一个线程块。这样将为每个线程块创建Mds和Nds数组，并且一个线程块的所有线程都可以访问相同的Mds和Nds。因为在一个块中的所有线程都必须能够访问由其相邻线程加载到Mds和Nds中的M和N元素，以便它们可以使用这些值来满足其输入需求。
第07行和第08行将threaddx和blockIdx值保存到名字较短的自动变量中，以使代码更简洁。回想一下标量变量将被放入寄存器中。它们的作用域属于每个线程。也就是说，运行时系统为每个线程都创建了tx、ty、bx和by的一个私有版本，它们将驻留在线程可访问的寄存器中。他们使用 threadIdx 和 blockIdx 值初始化，并在线程的生命周期内多次使用。一旦线程结束，这些变量的值不存在。
第 11 和 12 行分别确定线程要生成 P 元素的行索引和列索引。代码假定每个线程负责
计算一个 P 元素。如第 12 行所示，水平 (x) 位置或线程生成的 P 元素的列索引可以计算为 bx*TILE_WIDTH+tx。这是因为每个程序块在水平维度上覆盖了 P 的 TILE_WIDTH 元素。位于块 bx 中的线程之前会有 bx 个线程块，或 (bxTILE_WIDTH) 个线程；它们覆盖了 P 的 bx*TILE_WIDTH 个元素。因此，bx 和 tx 的线程应负责计算 x 索引为 bxTILE_WIDTH+tx 的 P 元素。在图 5.7 的例子中，程序块 1,0 的线程 0,1 计算的 P 元素的水平（x）索引为 02+1=1。该水平索引保存在线程的变量 Col 中，如下图所示。
类似地，线程要处理的P元素的垂直（y）位置或行索引计算为by*TILE_WIDTH+ty。在这个矩阵乘法的例子中，块1,0的线程0,1要计算的P元素的y索引是1*2+0=2。这个垂直索引保存在线程的变量Row中。如图5.10所示，每个线程计算Colth列和Rowth行的P元素。因此，块1,0的线程0,1要计算的P元素是P[2,1]。
上图的第16行标志着循环的开始，该循环遍历计算P元素的所有阶段。循环的每次迭代对应于之前图中的一个计算阶段。ph变量表示已经完成的点积阶段数。请记住，每个阶段使用一个M的块和一个N的块。因此，在每个阶段开始时，ph * TILE_WIDTH对M和N元素已经被前面的阶段处理过。
在每个阶段，上图中的第19和20行分别将适当的M和N元素加载到共享内存中。由于我们已经知道线程要处理的M的行和N的列，现在我们将注意力转向M的列索引和N的行索引。如下图所示，每个块有TILE_WIDTH²个线程，将协作加载TILE_WIDTH²个M元素和TILE_WIDTH²个N元素到共享内存中。因此，我们需要做的就是分配每个线程加载一个M元素和一个N元素。这可以方便地通过使用blockIdx和threadIdx来完成。请注意，要加载的M元素部分的起始列索引是ph*TILE_WIDTH。因此，一个简单的方法是让每个线程加载一个距离该起始点tx（值）位置的元素。同样，要加载的N元素部分的起始行索引也是ph*TILE_WIDTH。因此，每个线程加载一个距离该起始点ty（值）位置的元素。
这正是我们在第19和第20行中所看到的。在第19行，每个线程加载M[Row*Width + ph*TILE_WIDTH + tx]，其中线性索引是由行索引Row和列索引ph*TILE_WIDTH + tx形成的。由于Row的值是ty的线性函数，每个TILE_WIDTH²线程将把一个独特的M元素加载到共享内存中，因为每个线程都有一个独特的tx和ty组合。总的来说，这些线程将在上图加载M的一个暗方块子集。以类似的方式，在第20行中，每个线程使用线性索引(ph*TILE_WIDTH + ty)*Width + Col将适当的N元素加载到共享内存中。读者应使用前面图中的小示例来验证地址计算是否正确适用于各个线程。
__syncthreads() 在第21行确保所有线程在任何线程继续之前，已经完成将M和N的块加载到Mds和Nds中。回想一下CUDA编程04 - GPU计算架构和线程调度 中，用 __syncthreads() 可以使一个块中的所有线程相互等待，直到它们都到达同步点，然后才能继续。这一点很重要，因为线程使用的M和N元素可能会被其他线程加载。需要确保所有元素都正确加载到共享内存中，然后才能开始使用这些元素。第23行的循环则基于块元素执行点积的一个阶段。线程ty, tx的循环进程如上图所示，M和N元素的访问方向沿着标记为k的箭头，k是第23行的循环变量。请注意，这些元素将从Mds和Nds中访问，这些是存放这些M和N元素的共享内存数组。第26行的屏障 __syncthreads() 确保所有线程在继续到下一次迭代并从下一个块加载元素之前，已经完成了对共享内存中M和N元素的使用。因此，没有任何线程会过早加载元素，从而破坏其他线程的输入值。
在第21行和第26行中的两个__syncthreads()调用展示了并行程序员在协调线程时常常需要考虑的两种不同类型的数据依赖。第一种被称为读后写依赖，因为线程必须等待其他线程将数据写入适当的位置，然后才能尝试读取。第二种被称为写后读依赖，因为线程必须等待所有需要该数据的线程读取它，然后才能覆盖它。读后写和写后读依赖的其他名称分别是真依赖和假依赖。读后写依赖是真依赖，因为读取线程确实需要写入线程提供的数据，因此别无选择，只能等待。写后读依赖是假依赖，因为写入线程不需要读取线程的任何数据。依赖的产生是因为它们重用了相同的内存位置，如果使用不同的位置则不会存在这种依赖。
从第16行到第28行的循环嵌套展示了一种称为stripmining的技术，它将一个长时间运行的循环分解为多个阶段。每个阶段涉及一个内部循环，该循环执行原始循环的几个连续迭代。原始循环变成一个外部循环，其作用是迭代地调用内部循环，以便以原始顺序执行原始循环的所有迭代。通过在内部循环之前和之后添加栅栏同步，我们强制同一个块中的所有线程在每个阶段集中工作在同一输入数据部分。stripmining是创建数据并行程序中所需阶段的重要手段。
在点积的所有阶段完成后，执行退出外部循环。在第29行，所有线程使用从行和列计算出的线性索引来写入对应P中的元素。
分块算法的好处是显著的。在矩阵乘法中，全球内存访问减少了 TILE_WIDTH 的一个因子。使用 16 x 16 的平铺，可以将全球内存访问减少一个因子为 16。这将计算与全球内存访问的比率从 0.25 OP/B 提高到 4 OP/B。这一改进使得 CUDA 设备的内存带宽能够支持更高的计算速率。例如，在 A100 GPU 中，其全球内存带宽为 1555 GB/秒，这一改进使得设备能够达到 (1555 GB/秒) * (4 OP/B) = 6220 GFLOPS，这远高于未使用平铺的内核所达到的 389 GFLOPS。
尽管分块显著提高了吞吐量，但 6220 GFLOPS 仍然仅占设备峰值吞吐量 19,500 GFLOPS 的 32%。可以进一步优化代码以减少全球内存访问并提高吞吐量。我们将在本书后面看到一些这些优化，而其他高级优化将不予讨论。由于矩阵乘法在许多领域的重要性，已经有高度优化的库，如 cuBLAS 和 CUTLASS，这些库已经包含了许多这些高级优化。程序员可以使用这些库在他们的线性代数应用中立即实现接近峰值性能。
在提高矩阵乘法的吞吐量方面，特别是在应用程序方面，分块的有效性并不是唯一针对GPU的。长期以来，应用分块技术以提高CPU性能的历史悠久，其方法是确保在特定时间窗口内被CPU线程重用的数据能够在缓存中找到。一个关键的区别是，CPU上的分块技术依赖于CPU缓存隐式地将重用数据保留在芯片上，而GPU上的分块技术则显式地使用共享内存将数据保留在芯片上。原因在于，CPU核心通常同时运行一到两个线程，因此一个线程可以依赖缓存保持最近使用的数据。相比之下，GPU的流处理器（SM）同时运行多个线程，以便能够隐藏延迟。这些线程可能会争夺缓存槽，这使得GPU缓存的可靠性降低，因此需要使用共享内存来存储重要的重用数据。
尽管分块矩阵乘法核函数的性能提升令人印象深刻，但它确实做出了一些简化假设。首先，假设矩阵的宽度是线程块宽度的倍数。这使得内核无法正确处理任意宽度的矩阵。第二个假设是矩阵是方阵。这在实际中并不总是成立。在下一节中，我们将介绍一个带有边界检查的核函数，以消除这些假设。
七：边界处理
我们现在扩展分块矩阵乘法的核函数来处理任意宽度的矩阵。扩展必须允许核函数正确处理宽度不是分块宽度的倍数的矩阵。让我们使用 3x3 的M、N和P矩阵。修改后的示例如下图所示，请注意，矩阵的宽度为3，而不是分块宽度（2）的倍数。下图显示出了在block[0,0]的内存访问情况，我们看到Thread[0,1]和Thread[1,1]将尝试加载M中不存在的元素。类似地，我们看到线程[1，0]和线程[1,1]将尝试访问N中不存在的元素。
访问不存在的元素在两个方面都是有问题的。在访问超出行末尾的不存在元素的情况下（上图中的thread[0,1]和thread[1,1]对M的访问），这些访问将会访问到错误的元素。在我们的示例中，线程将尝试访问M[0,3]和M[1,3]，这些元素并不存在。那么这些内存加载会发生什么呢？为了回答这个问题，我们需要回到二维矩阵的线性化布局。在线性化布局中，M[0,2]之后的元素是M[1,0]。尽管thread[0,1]正在尝试访问M[0,3]，但最终将得到M[1,0]。在随后的内积计算中使用这个值显然会破坏输出值。
在访问超出列末尾的元素时也会出现类似的问题（上图中的thread[1,0]和thread[1,1]对N的访问）。这些访问是对数组分配区域之外的内存位置的访问。在某些系统中，它们将返回来自其他数据结构的随机值。在其他系统中，这些访问将被拒绝，导致程序中止。无论哪种情况，这种访问的结果都是不可取的。
根据我们到目前的讨论，问题访问似乎仅在线程执行的最后阶段出现。这表明我们可以通过在分块核函数执行的最后阶段采取特殊措施来处理它。不幸的是，这并不正确。有问题的访问可以在所有阶段出现。下图显示了在阶段0中block[1,1]的内存访问情况。我们看到thread[1,0]和thread[1,1]试图访问不存在的M元素M[3,0]和M[3,1]，而thread[0,1]和thread[1,1]试图访问不存在的N[0,3]和N[1,3]。
请注意，这些问题访问不能简单的通过排除线程的方式解决。例如，block[1,1] 中的 thread[1,0] 不计算任何有效的 P 中元素。然而，它需要在阶段 0 加载 M[2,1]，以供 block[1,1] 中的其他线程使用。此外，请注意，一些计算有效 P 中元素的线程将尝试访问不存在的 M 或 N 元素。例如，正如我们在上图中看到的，block [0,0] 的 thread[0,1] 计算了一个有效的 P 元素 P[0,1]。然而，它在阶段 1 尝试访问一个不存在的 M[0,3]。这两个事实表明，我们需要使用不同的边界条件测试, 来加载 M 、加载 N 以及计算/存储 P 中元素。一个经验法则是，每个内存访问都需要有一个相应的检查，以确保访问中使用的索引在被访问数组的范围内。
让我们从加载分块矩阵的边界测试条件开始。当一个线程要加载一个分块矩阵元素时，它应该测试它试图加载的输入元素是否是有效元素。这可以通过检查 y 和 x 索引轻松完成。例如，在上图代码的第 19 行，线性化索引是由行的 y 索引和 ph*TILE_WIDTH + tx 的 x 索引计算出来的。边界条件测试将是两个索引都小于 Width: Row < Width && (ph*TILE_WIDTH+tx) < Width。如果条件为真，线程应该继续加载 M 元素。读者应该验证加载 N 元素的条件测试是 (ph*TILE_WIDTH + ty) < Width && Col < Width。
如果条件为假，线程不应该加载该元素。问题是共享内存位置的应该放入什么。答案是 0.0，这是一个在内积计算中使用时不会造成任何副作用的值。如果任何线程在其内积计算中使用这个 0.0 值，内积值将不会发生任何变化。
最后，线程只有在负责计算有效的 P 元素时才应存储其最终的内积值。测试此条件的方式是 (Row < Width) && (Col < Width)。附加边界条件检查的内核代码如下图所示。
通过边界条件检查，块矩阵乘法核函数离成为通用矩阵乘法内核仅一步之遥。一般来说，矩阵乘法是定义在长方形矩阵上的：一个 j x k 的 M 矩阵与一个 k x l 的 N 矩阵相乘，结果是一个 j x l 的 P 矩阵。到目前为止，我们的内核只能处理方阵。
幸运的是，将我们的内核进一步扩展为通用矩阵乘法内核是相当简单的。我们需要做一些简单的更改。首先，Width 参数被三个无符号整数参数替代：j、k、l。将 Width 用于指代 M 的高度或 P 的高度的地方，替换为 j。将 Width 用于指代 M 的宽度或 N 的高度的地方，替换为 k。将 Width 用于指代 N 的宽度或 P 的宽度的地方，替换为 l。对核函数进行这些更改请读者自行练习。
八：内存使用对占用率的影响（occupancy）
回顾一下在 CUDA编程04 - GPU计算架构和线程调度中，我们讨论了最大化线程在SM上的占用率以容忍长延迟操作的重要性。核函数的内存使用在占用率调优中起着重要作用。虽然CUDA寄存器和共享内存在减少对全局内存访问次数方面非常有效，但必须小心这些SM的内存的容量限制。每个CUDA设备提供的资源是有限的，这限制了在给定应用程序中可以同时驻留在SM中的线程数量。一般来说，每个线程所需的资源越多，可以驻留在每个SM中的线程数量就越少。
我们在 CUDA编程04 - GPU计算架构和线程调度 中看到，寄存器使用可能是占用率的限制因素。共享内存使用也可能限制可以分配给每个SM的线程数量。例如，A100 GPU可以配置为每个SM最多有164 KB的共享内存，并支持每个SM最多2048个线程。因此，为了使用所有2048个线程槽，线程块的平均使用量不应超过（164 KB）/（2048线程）= 82 B/线程。在平铺矩阵乘法示例中，每个块有TILE_WIDTH²个线程，并使用TILE_WIDTH²*4B的共享内存用于Mds，以及TILE_WIDTH²*4B的共享内存用于Nds。因此，线程块的共享内存平均使用量为（TILE_WIDTH²*4B + TILE_WIDTH²*4B）/（TILE_WIDTH²线程）= 8 B/线程的共享内存。因此分块内存核函数的占用率不受共享内存限制。
然而，考虑一个核函数，其线程块使用32 KB的共享内存，每个线程块有256个线程。在这种情况下，内核平均使用(32 KB)/(256线程)=132 B/线程的共享内存。由于这样的共享内存使用，核函数无法实现最大占用率。每个SM最多只能容纳(164 KB)/(132 B/线程)=1272个线程。因此，该核的函数最大可实现占用率将是(1272分配线程)/(2048最大线程)=62%。
请注意，每个SM中的共享内存大小可能因设备而异。每一代或型号的设备在每个SM中可能具有不同数量的共享内存。通常希望内核能够根据硬件中可用的共享内存量使用不同数量的共享内存。也就是说，我们可能希望主机代码动态确定共享内存的大小，并调整内核使用的共享内存量。这可以通过调用cudaGetDeviceProperties函数来完成。假设变量&devProp被传递给该函数。在这种情况下，字段给出了每个SM中可用的共享内存量。程序员可以据此确定每个块应使用的共享内存量。
不幸的是，我们例子的核函数不支持主机代码对共享内存使用的任何动态调整。在上图中使用的声明将其共享内存使用的大小硬编码为编译时常量。
也就是说，Mds和Nds的大小被设置为TILE_WIDTH * TILE_WIDTH个元素，无论编译时TILE_WIDTH的值设置何值。代码是一个常量值：
Mds和Nds都将具有256个元素。如果我们想要更改Mds和Nds的大小，我们需要更改TILE_WIDTH的值并重新编译代码。内核在运行时无法轻松调整其共享内存使用，而无需重新编译。
我们可以通过在共享内存声明前添加C extern关键字并在声明中省略数组大小，使用CUDA中的另一种声明风格来启用这种调整。基于这种风格，Mds和Nds的声明需要合并为一个动态分配的数组。
由于只有一个动态数组，我们还需要手动定义数组的 Mds 部分开始的位置和 Nds 部分开始的位置。请注意，动态数组是一维的。我们需要通过基于垂直和水平索引的线性化索引来访问它。
在运行时，即当我们调用核函数时，我们可以根据设备查询结果动态配置每个块使用的共享内存大小，并将其作为第三个配置参数提供给核函数调用。例如，修该后的核函数可以通过以下语句启动
其中 size_t 是内置类型，用于声明变量以保存动态分配数据结构的大小信息的。大小以字节数表示。在我们的矩阵乘法示例中，对于一个 16 x 16 的瓦片，我们有一个大小为 2 x 16 x 16 x 4=2048 字节，以容纳 Mds 和 Nds。
在下图中，我们展示了修改后的内核代码，以使用动态大小的共享内存来存储Mds和Nds数组。将数组每个部分的大小作为参数传递给核函数。在这个示例中，我们添加了两个参数：第一个参数是Mds部分的大小，第二个参数是Nds部分的大小，单位为字节。请注意，在上述主机代码中，我们将size/2作为这些参数的值传递，值为1024字节。通过第06行和第07行的赋值，核函数代码的其余部分可以将Mds和Nds用作数组的基址，并使用线性索引访问Mds和Nds元素。例如，使用Mds[ty][tx]时，可以改为使用Mds[ty*TILE_WIDTH+tx]。
九：总结
总之，现代处理器中程序的执行速度可能会受到内存速度的严重限制。为了实现CUDA设备的执行吞吐量的良好利用，需要在核函数代码中努力提高计算与全局内存访问的比率。如果比率较低，则内核会受到内存的限制。也就是说，其执行速度受到从内存访问操作数的速率的限制。
CUDA提供了对寄存器、共享内存和常量内存的访问。这些内存比全局内存小得多，但可以以更高的速度访问。有效使用这些内存需要重新设计算法。我们以矩阵乘法为例，说明了切块（tiling），这是一种流行的策略，用于增强数据访问的局部性并有效利用共享内存。在并行编程中，切块使用屏障同步强制多个线程在执行的每个阶段共同关注输入数据的一个子集，以便将该子集数据放入这些特殊内存类型中，从而实现更高的访问速度。
然而，CUDA程序员需要意识到这些特殊类型内存的大小有限。它们的容量依赖于具体的实现。一旦超出其容量，就会限制每个SM中可以同时执行的线程数量，并可能对GPU的计算吞吐量以及其容忍延迟的能力产生负面影响。在开发应用程序时，能够考虑硬件限制是并行编程的一个关键方面。
尽管我们在CUDA C编程的背景下引入了分块算法，但这是一种在几乎所有类型的并行计算系统中实现高性能的有效策略。原因在于，应用程序必须表现出数据访问的局部性，以有效利用这些系统中的高速内存。例如，在多核CPU系统中，数据局部性使得应用程序能够有效使用片上数据缓存，从而减少内存访问延迟并实现高性能。这些片上数据缓存的大小也有限，并且需要计算表现出局部性。因此，读者在使用其他编程模型为其他类型的并行计算系统开发并行应用程序时，也会发现分块算法非常有用。
我们在本篇文章的目标是介绍局部性、分块和不同的CUDA内存类型的概念。我们介绍了使用共享内存的分块矩阵乘法核函数。我们进一步研究了边界测试条件的必要性，以便在应用分块技术时允许任意数据维度。我们还简要讨论了动态大小的共享内存分配的使用，以便内核可以根据硬件能力调整每个块使用的共享内存的大小。我们没有讨论在分块中使用寄存器。我们将在后面文章中讨论在分块算法中使用寄存器的情况。