GPU的线程模型和内存模型 遇见C++ AMP：在GPU上做并行计算 遇见C++ AMP：GPU的线程模型和内存模型
Written by Allen Lee
I see all the young believers, your target audience. I see all the old deceivers; we all just sing their song.
– Marilyn Manson, Target Audience (Narcissus Narcosis)
从CPU到GPU
在《遇见C++ PPL：C++的并行和异步》里，我们介绍了如何使用C++ PPL在CPU上做并行计算，这次，我们会把舞台换成GPU，介绍如何使用C++ AMP在上面做并行计算。
为什么选择在GPU上做并行计算呢？现在的多核CPU一般都是双核或四核的，如果把超线程技术考虑进来，可以把它们看作四个或八个逻辑核，但现在的GPU动则就上百个核，比如中端的NVIDIA GTX 560 SE就有288个核，*的NVIDIA GTX 690更有多达3072个核，这些超多核（many-core）GPU非常适合大规模并行计算。
接下来，我们将会在《遇见C++ PPL：C++的并行和异步》的基础上，对并行计算正弦值的代码进行一番改造，使之可以在GPU上运行。如果你没读过那篇文章，我建议你先去读一读它的第一节。此外，本文也假设你对C++ Lambda有所了解，否则，我建议你先去读一读《遇见C++ Lambda》。
并行计算正弦值
首先，包含/引用相关的头文件/命名空间，如代码1所示。amp.h是C++ AMP的头文件，包含了相关的函数和类，它们位于concurrency命名空间之内。amp_math.h包含了常用的数学函数，如sin函数，concurrency::fast_math命名空间里的函数只支持单精度浮点数，而concurrency::precise_math命名空间里的函数则对单精度浮点数和双精度浮点数均提供支持。
代码 1
把浮点数的类型从double改成float，如代码2所示，这样做是因为并非所有GPU都支持双精度浮点数的运算。另外，std和concurrency两个命名空间都有一个array类，为了消除歧义，我们需要在array前面加上"std::"前缀，以便告知编译器我们使用的是STL的array类。
代码 2
接着，创建一个array_view对象，把前面创建的array对象包装起来，如代码3所示。array_view对象只是一个包装器，本身不能包含任何数据，必须和真正的容器搭配使用，如C风格的数组、STL的array对象或vector对象。当我们创建array_view对象时，需要通过类型参数指定array_view对象里的元素的类型以及它的维度，并通过构造函数的参数指定对应维度的长度以及包含实际数据的容器。
代码 3
代码3创建了一个一维的array_view对象，这个维度的长度和前面的array对象的长度一样，这个包装看起来有点多余，为什么要这样做？这是因为在GPU上运行的代码无法直接访问系统内存里的数据，需要array_view对象出来充当一个桥梁的角色，使得在GPU上运行的代码可以通过它间接访问系统内存里的数据。事实上，在GPU上运行的代码访问的并非系统内存里的数据，而是复制到显存的副本，而负责把这些数据从系统内存复制到显存的正是array_view对象，这个过程是自动的，无需我们干预。
有了前面这些准备，我们就可以着手编写在GPU上运行的代码了，如代码4所示。parallel_for_each函数可以看作C++ AMP的入口点，我们通过extent对象告诉它创建多少个GPU线程，通过Lambda告诉它这些GPU线程运行什么代码，我们通常把这个代码称作Kernel。
代码 4
我们希望每个GPU线程可以完成和结果集里的某个元素对应的一组操作，比如说，我们需要计算10个浮点数的正弦值，那么，我们希望创建10个GPU线程，每个线程依次完成读取浮点数、计算正弦值和保存正弦值三个操作。但是，每个GPU线程运行的代码都是一样的，如何区分不同的GPU线程，并定位需要处理的数据呢？
这个时候就轮到index对象出场了，我们的array_view对象是一维的，因此index对象的类型是index<1>，这个维度的长度是10，因此将会产生从0到9的10个index对象，每个GPU线程对应其中一个index对象。这个index对象将会通过Lambda的参数传给我们，而我们将会在Kernel里通过这个index对象找到当前GPU线程需要处理的数据。
既然Lambda的参数只传递index对象，那Kernel又是如何与外界交换数据的呢？我们可以通过闭包捕获当前上下文的变量，这使我们可以灵活地操作多个数据源和结果集，因此没有必要提供返回值。从这个角度来看，C++ AMP的parallel_for_each函数在用法上类似于C++ PPL的parallel_for函数，如代码5所示，我们传给前者的extent对象代替了我们传给后者的起止索引值。
代码 5
那么，Kernel右边的restrict(amp)修饰符又是怎么一回事呢？Kernel最终是在GPU上运行的，不管以什么样的形式，restrict(amp)修饰符正是用来告诉编译器这点的。当编译器看到restrict(amp)修饰符时，它会检查Kernel是否使用了不支持的语言特性，如果有，编译过程中止，并列出错误，否则，Kernel会被编译成HLSL，并交给DirectCompute运行。Kernel可以调用其他函数，但这些函数必须添加restrict(amp)修饰符，比如代码4的sin函数。
计算完毕之后，我们可以通过一个for循环输出array_view对象的数据，如代码6所示。当我们在CPU上首次通过索引器访问array_view对象时，它会把数据从显存复制回系统内存，这个过程是自动的，无需我们干预。
代码 6
哇，不知不觉已经讲了这么多，其实，使用C++ AMP一般只涉及到以下三步：
- 创建array_view对象。
- 调用parallel_for_each函数。
- 通过array_view对象访问计算结果。
其他的事情，如显存的分配和释放、GPU线程的规划和管理，C++ AMP会帮我们处理的。
并行计算矩阵之和
上一节我们通过一个简单的示例了解C++ AMP的使用步骤，接下来我们将会通过另一个示例深入了解array_view、extent和index在二维场景里的用法。
假设我们现在要计算两个100 x 100的矩阵之和，首先定义矩阵的行和列，然后通过create_matrix函数创建两个vector对象，接着创建一个vector对象用于存放矩阵之和，如代码7所示。
代码 7
create_matrix函数的实现很简单，它接受矩阵的总容量（行和列之积）作为参数，然后创建并返回一个包含100以内的随机数的vector对象，如代码8所示。
代码 8
值得提醒的是，当create_matrix函数执行"return matrix;"时，会把vector对象拷贝到一个临时对象，并把这个临时对象返回给调用方，而原来的vector对象则会因为超出作用域而自动销毁，但我们可以通过编译器的Named Return Value Optimization对此进行优化，因此不必担心按值返回会带来性能问题。
虽然我们通过行和列等二维概念定义矩阵，但它的实现是通过vector对象模拟的，因此在使用的时候我们需要做一下索引变换，矩阵的第m行第n列元素对应的vector对象的索引是m * columns + n（m、n均从0开始计算）。假设我们要用vector对象模拟一个3 x 3的矩阵，如图1所示，那么，要访问矩阵的第2行第0列元素，应该使用索引6（2 * 3 + 0）访问vector对象。
图 1
接下来，我们需要创建三个array_view对象，分别包装前面创建的三个vector对象，创建的时候先指定行的大小，再指定列的大小，如代码9所示。
代码 9
因为我们创建的是二维的array_view对象，所以我们可以直接使用二维索引访问矩阵的元素，而不必像前面那样计算对应的索引。还是以3 x 3的矩阵为例，如图2所示，vector对象会被分成三段，每段包含三个元素，第一段对应array_view对象的第一行，第二段对应第二行，如此类推。如果我们想访问矩阵的第2行第0列的元素，可以直接使用索引 (2, 0) 访问array_view对象，这个索引对应vector对象的索引6。
图 2
考虑到第一、二个array_view对象的数据流动方向是从系统内存到显存，我们可以把它们的第一个类型参数改为const int，如代码10所示，表示它们在Kernel里是只读的，不会对它包装的vector对象产生任何影响。至于第三个array_view对象，由于它只是用来输出计算结果，我们可以在调用parallel_for_each函数之前调用array_view对象的discard_data成员函数，表明我们对它包装的vector对象的数据不感兴趣，不必把它们从系统内存复制到显存。
代码 10
有了这些准备，我们就可以着手编写Kernel了，如代码11所示。我们把第三个array_view对象的extent传给parallel_for_each函数，由于这个矩阵是100 x 100的，parallel_for_each函数会创建10,000个GPU线程，每个GPU线程计算这个矩阵的一个元素。由于我们访问的array_view对象是二维的，索引的类型也要改为相应的index<2>。
代码 11
看到这里，你可能会问，GPU真能创建这么多个线程吗？这取决于具体的GPU，比如说，NVIDIA GTX 690有16个多处理器（Kepler架构，每个多处理器有192个CUDA核），每个多处理器的最大线程数是2048，因此可以同时容纳最多32,768个线程；而NVIDIA GTX 560 SE拥有9个多处理器（Fermi架构，每个多处理器有32个CUDA核），每个多处理器的最大线程数是1536，因此可以同时容纳最多13,824个线程。
计算完毕之后，我们可以在CPU上通过索引器访问计算结果，代码12向控制台输出结果矩阵的第14行12列元素。
代码 12
async + continuation
掌握了C++ AMP的基本用法之后，我们很自然就想知道parallel_for_each函数会否阻塞当前CPU线程。parallel_for_each函数本身是同步的，它负责发起Kernel的运行，但不会等到Kernel的运行结束才返回。以代码13为例，当parallel_for_each函数返回时，即使Kernel的运行还没结束，checkpoint 1位置的代码也会照常运行，从这个角度来看，parallel_for_each函数是异步的。但是，当我们通过array_view对象访问计算结果时，如果Kernel的运行还没结束，checkpoint 2位置的代码会卡住，直到Kernel的运行结束，array_view对象把数据从显存复制到系统内存为止。
代码 13
既然Kernel的运行是异步的，我们很自然就会希望C++ AMP能够提供类似C++ PPL的continuation。幸运的是，array_view对象提供一个synchronize_async成员函数，它返回一个concurrency::completion_future对象，我们可以通过这个对象的then成员函数实现continuation，如代码14所示。事实上，这个then成员函数就是通过C++ PPL的task对象实现的。
代码 14
你可能会问的问题
1. 开发C++ AMP程序需要什么条件？
你需要Visual Studio 2012以及一块支持DirectX 11的显卡，Visual C++ 2012 Express应该也可以，如果你想做GPU调试，你还需要Windows 8操作系统。运行C++ AMP程序需要Windows 7/Windows 8以及一块支持DirectX 11的显卡，部署的时候需要把C++ AMP的运行时（vcamp110.dll）放在程序可以找到的目录里，或者在目标机器上安装Visual C++ 2012 Redistributable Package。
2. C++ AMP是否支持其他语言？
C++ AMP只能在C++里使用，其他语言可以通过相关机制间接调用你的C++ AMP代码：
- How to use C++ AMP from C#
- How to use C++ AMP from C# using WinRT
- How to use C++ AMP from C++ CLR app
- Using C++ AMP code in a C++ CLR project
3. C++ AMP是否支持其他平台？
目前C++ AMP只支持Windows平台，不过，微软发布了C++ AMP开放标准，支持任何人在任何平台上实现它。如果你希望在其他平台上利用GPU做并行计算，你可以考虑其他技术，比如NVIDIA的CUDA（只支持NVIDIA的显卡），或者OpenCL，它们都支持多个平台。
4. 能否推荐一些C++ AMP的学习资料？
目前还没有C++ AMP的书，Kate Gregory和Ade Miller正在写一本关于C++ AMP的书，希望很快能够看到它。下面推荐一些在线学习资料：
- C++ AMP open specification
- Parallel Programming in Native Code (team blog)
- C++ AMP (C++ Accelerated Massive Parallelism)
- C++ AMP Videos
*声明：本文已经首发于InfoQ中文站，版权所有，《遇见C++ AMP：在GPU上做并行计算》，如需转载，请务必附带本声明，谢谢。
遇见C++ AMP：GPU的线程模型和内存模型
Written by Allen Lee
I don't care where the enemies are / Can't be stopped / All I know / Go hard
– Linkin Park, Lost In The Echo
C++ AMP、CUDA和OpenCL，选择哪个？
在《遇见C++ AMP：在GPU上做并行计算》发布之后，我曾被多次问及为何选择C++ AMP，以及它与CUDA、OpenCL等相比有何优势，看来有必要在进入正题之前就这个问题发表一下看法了。
在众多可以影响决策的因素之中，平台种类的支持和GPU种类的支持是两个非常重要的因素，它们联合起来足以直接否决某些选择。如果我们把这两个因素看作两个维度，可以把平面分成四个象限，C++ AMP、CUDA和OpenCL分别位于第二象限、第四象限和第一象限，如图1所示。如果你想通吃所有平台和所有GPU，OpenCL是目前唯一的选择，当然，你也需要为此承担相当的复杂性。CUDA是一个有趣的选择，紧贴最新的硬件技术、数量可观的行业应用和类库支持使之成为一个无法忽视的选择，但是，它只能用于NVIDIA的GPU极大地限制了它在商业应用上的采用，我想你不会为了运行我的应用特意把显卡换成NVIDIA的。C++ AMP的情况刚好相反，它适用于各种支持DirectX 11的GPU，但只能在Windows上运行。
图 1
这些技术都有自己的特点和位置，你应该根据项目的具体情况选择合适的解决方案。如果你正在从事的工作需要进行大量计算，你想尽可能利用硬件特性对算法进行优化，而你的机器刚好有一块NVIDIA的显卡，并且你不需要在其他机器上重复执行这些计算，那么CUDA将是你的不二之选。尽管NVIDIA已经开源CUDA编译器，并且欢迎其他厂商通过CUDA编译器SDK添加新的语言/处理器，但AMD不太可能会为它提供在AMD的GPU上运行的扩展，毕竟它也有自己的基于OpenCL的AMD APP技术。如果你正在从事Windows应用程序的开发工作，熟悉C++和Visual Studio，并且希望借助GPU进一步提升应用程序的性能，那么C++ AMP将是你的不二之选。尽管微软已经开放C++ AMP规范，Intel的Dillon Sharlet也通过Shevlin Park项目验证了在Clang/LLVM上使用OpenCL实现C++ AMP是可行的，但这不是一个产品级别的商用编译器，Intel也没有宣布任何发布计划。如果你确实需要同时兼容Windows、Mac OS X和Linux等多个操作系统，并且需要同时支持NVIDIA和AMD的GPU，那么OpenCL将是你的不二之选。
GPU线程的执行
在《遇见C++ AMP：在GPU上做并行计算》里，我们通过extent对象告诉parallel_for_each函数创建多少个GPU线程，那么，这些GPU线程又是如何组织、分配和执行的呢？
首先，我们创建的GPU线程会被分组，分组的规格并不固定，但必须满足两个条件：对应的维度必须能被整除，分组的大小不能超过1024。假设我们的GPU线程是一维的，共8个，如图2所示，则可以选择每2个GPU线程为1组或者每4个GPU线程为1组，但不能选择每3个GPU线程为1组，因为剩下的2个GPU线程不足1组。
图 2
假设我们创建的GPU线程是二维的，3 x 4，共12个，如图3所示，则可以选择3 x 1或者3 x 2作为分组的规格，但不能选择2 x 2作为分组的规格，因为剩下的4个GPU线程虽然满足分组的大小，但不满足分组的形状。每个分组必须完全相同，包括大小和形状。
图 3
为了便于解释，我们的GPU线程只有寥寥数个，但真实案例的GPU线程往往是几十万甚至几百万个，这个时候，分组的规格会有大量选择，我们必须仔细判断它们是否满足条件。假设我们的GPU线程是640 x 480，那么16 x 48、32 x 16和32 x 32都可以选择，它们分别产生40 x 10、20 x 30和20 x 15个分组，但32 x 48不能选择，因为它的大小已经超过1024了。
接着，这些分组会被分配到GPU的流多处理器（streaming multiprocessor），每个流多处理器根据资源的使用情况可能分得一组或多组GPU线程。在执行的过程中，同一组的GPU线程可以同步，不同组的GPU线程无法同步。你可能会觉得这种有限同步的做法会极大地限制GPU的作为，但正因为组与组之间是相互独立的，GPU才能随意决定这些分组的执行顺序。这有什么好处呢？假设低端的GPU每次只能同时执行2个分组，那么执行8个分组需要4个执行周期，假设高端的GPU每次可以同时执行4个分组，执行8个分组只需2个执行周期，如图4所示，这意味着我们写出来的程序具备可伸缩性，能够自动适应GPU的计算资源。
图 4
说了这么多，是时候看看代码了。parallel_for_each函数有两种模式，一种是简单模式，我们通过extent对象告诉它创建多少GPU线程，C++ AMP负责对GPU线程进行分组，另一种是分组模式，我们通过tiled_extent对象告诉它创建多少GPU线程以及如何进行分组。创建tiled_extent对象非常简单，只需在现有的extent对象上调用tile方法，并告知分组的规格就行了，如代码1所示。值得提醒的是，分组的规格是通过模板参数告诉tile方法的，这意味着分组的规格必须在编译时确定下来，C++ AMP目前无法做到运行时动态分组。
代码 1
既然C++ AMP不支持运行时动态分组，肯定会为简单模式预先定义一些分组的规格，那么C++ AMP又是如何确保它们能被整除？假设我们创建的GPU线程是一维的，共10000个，C++ AMP会选择每256个GPU线程为1组，把前面9984个GPU线程分成39个分组，然后补充240个GPU线程和剩下的16个GPU线程凑够1组，执行的时候会通过边界测试确保只有前10000个GPU线程执行我们的代码。对于二维和三维的情况，C++ AMP也会采取这种补充GPU线程的策略，只是分组的规格不同，必要时还会重新排列GPU线程，以便分组能够顺利完成。需要说明的是，简单模式背后采取的策略属于实现细节，在这里提及是为了满足部分读者的好奇心，你的算法不该对它有所依赖。
共享内存的访问
既然简单模式可以自动分组，为何还要大费周章使用分组模式？为了回答这个问题，我们先要了解一下GPU的内存模型。在Kernel里，我们可以访问全局内存、共享内存和寄存器，如图5所示。当我们通过array_view对象把数据从主机内存复制到显卡内存时，这些数据会被保存在全局内存，直到应用程序退出，所有GPU线程都能访问全局内存，不过访问速度很慢，大概需要1000个GPU时钟周期，大量的GPU线程反复执行这种高延迟的操作将会导致GPU计算资源的闲置，从而降低整体的计算性能。
图 5
为了避免反复从全局内存访问相同的数据，我们可以把这些数据缓存到寄存器或者共享内存，因为它们集成在GPU芯片里，所以访问速度很快。当我们在Kernel里声明一个基本类型的变量时，它的数据会被保存在寄存器，直到GPU线程执行完毕，每个GPU线程只能访问自己的寄存器，寄存器的容量非常小，不过访问速度非常快，只需1个GPU时钟周期。当我们在Kernel里通过tile_static关键字声明一个变量时，它的数据会被保存在共享内存（也叫tile_static内存），直到分组里的所有GPU线程都执行完毕，同一组的GPU线程都能访问相同的共享内存，共享内存的容量很小，不过访问速度很快，大概需要10个GPU时钟周期。tile_static关键字只能在分组模式里使用，因此，如果我们想使用共享内存，就必须使用分组模式。
如果数据只在单个GPU线程里反复使用，可以考虑把数据缓存到寄存器。如果数据会在多个GPU线程里反复使用，可以考虑把数据缓存到共享内存。共享内存的缓存策略是对全局内存的数据进行分组，然后把这些分组从全局内存复制到共享内存。假设我们需要缓存4 x 4的数据，可以选择2 x 2作为分组的规格把数据分成4组，如图6所示。以右上角的分组为例，我们需要4个GPU线程分别把这4个数据从全局内存复制到共享内存。复制的过程涉及两种不同的索引，一种是相对于所有数据的全局索引，用于从全局内存访问数据，另一种是相对于单个分组的本地索引，用于从共享内存访问数据，比如说，全局索引(1, 2)对应本地索引(1, 0)。
图 6
在分组模式里，我们可以通过tiled_index对象访问索引信息，它的global属性返回全局索引，local属性返回本地索引，tile属性返回分组索引，它是分组作为一个整体相对于其他分组的索引，tile_origin属性返回分组原点的全局索引，它是分组里的(0, 0)位置上的元素的全局索引。还是以右上角的分组为例，(1, 2)位置的global属性的值是(1, 2)，local属性的值是(1, 0)，tile属性的值是(0, 1)，tile_origin属性的值是(0, 2)。tiled_index对象将会通过Lambda的参数传给我们，我们将会在Kernel里通过它的属性访问全局内存和共享内存。
说了这么多，是时候看看代码了。正如extent对象搭配index对象用于简单模式，tiled_extent对象搭配tiled_index对象用于分组模式，使用的时候，两者的模板参数必须完全匹配，如代码2所示。parallel_for_each函数将会创建16个GPU线程，每4个GPU线程为1组，同一组的GPU线程共享一个2 x 2的数组变量，每个元素由一个GPU线程负责复制，每个GPU线程通过tiled_index对象的global属性获知从全局内存的哪个位置读取数据，通过local属性获知向共享内存的哪个位置写入数据。
代码 2
因为缓存的数据会在多个GPU线程里使用，所以每个GPU线程必须等待其他GPU线程缓存完毕才能继续执行后面的代码，否则，一些GPU线程还没开始缓存数据，另一些GPU线程就开始使用数据了，这样计算出来的结果肯定是错的。为了避免这种情况的发生，我们需要在代码2后面加上一句idx.barrier.wait();，加上之后的效果就像设了一道闸门，如图7所示，它把整个代码分成两个阶段，第一阶段缓存数据，第二阶段计算结果，缓存完毕的GPU线程会在闸门前面等待，当所有GPU线程都缓存完毕时，就会打开闸门让它们进入第二阶段。
图 7
总的来说，使用分组模式是为了借助共享内存减少全局内存的访问，缓存的过程已经包含了一次全局内存的访问，因此，如果我们的算法只需访问全局内存一次，比如《遇见C++ AMP：在GPU上做并行计算》的"并行计算矩阵之和"，那么缓存数据不会带来任何改善，反而增加了代码的复杂性。
并行计算矩阵之积
矩阵的乘法需要反复访问相同的元素，非常适合用来演示分组模式。接下来，我们将会分别使用简单模式和分组模式实现矩阵的乘法，然后通过对比了解这两种实现的区别。
设矩阵
求AB。设C = AB，根据定义，，其中，
。你可以把这个公式想象成矩阵A的第i行和矩阵B的第j列两个数组对应位置的元素相乘，然后相加。
如何把这些数学描述翻译成代码呢？第一步，定义A、B和C三个矩阵，如代码3所示，iota函数可以在指定的起止位置之间填充连续的数字，正好满足这里的需求。
代码 3
第二步，计算矩阵C的元素，如代码4所示，整个Kernel就是计算的求和公式， 因为每个元素的计算都是独立的，所以非常适合并行执行。
代码 4
在执行代码4的时候，parallel_for_each函数将会创建36个GPU线程，每个GPU线程计算矩阵C的一个元素，因为这36个GPU线程会同时执行，所以计算矩阵C的时间就是计算一个元素的时间。这听起来已经很好，还能更好吗？仔细想想，计算需要访问矩阵A的第i行一次，那么，计算矩阵C的第i行将会访问矩阵A的第i行M次，M是矩阵C的列数，在这里是6；同理，计算矩阵C的第j列将会访问矩阵B的第j列M次，M是矩阵C的行数，在这里也是6。因为A、B和C三个矩阵的数据是保存在全局内存的，所以优化的关键就是减少全局内存的访问。
根据上一节的讨论，我们将会使用分组模式，并把需要反复访问的数据从全局内存缓存到共享内存，那么，使用分组模式会对性能带来多少改善，又对算法造成多少影响呢，这正是我们接下来需要探讨的。
第一步，选择2 x 2作为分块的规格对A、B两个矩阵进行分块处理
分块矩阵的乘法和普通矩阵的乘法是一样的，设C = AB，根据定义，分块矩阵，其中，
。
第二步，把parallel_for_each函数改成分组模式，如代码5所示。T是子块的边长，W是分块矩阵A的列数，也是分块矩阵B的行数。
代码 5
第三步，分别缓存和
，如代码6所示。因为它们都是2 x 2的矩阵，所以缓存它们的工作需要4个GPU线程协同完成。正确缓存的关键在于弄清每个GPU线程负责全局内存和共享内存的哪些位置，共享内存的位置可以通过tiled_index对象的local属性获知，而全局内存的位置则需要换算一下，因为i和j是针对矩阵C而不是矩阵A和矩阵B的。每个GPU线程只是分别从矩阵A和矩阵B缓存一个元素，根据定义，从矩阵A缓存的元素必定位于第i行，而从矩阵B缓存的元素必定位于第j列。当我们缓存
时，子块位于分块矩阵A的左上角，tiled_index对象的local属性和global属性指向相同的列，因此，目标元素位于矩阵A的第h列，当我们缓存
时，我们已经从左到右跨过了w个子块，因此，目标元素位于矩阵A的第h + w * T列。同理，当我们缓存
时，子块位于分块矩阵B的左上角，目标元素位于矩阵B的第g行，当我们缓存
时，我们已经从上到下跨过了w个子块，因此，目标元素位于矩阵B的第g + w * T行。4个GPU线程都缓存完毕就会进入第二阶段。
代码 6
第四步，计算，如代码7所示。这是两个2 x 2的普通矩阵相乘，需要4个GPU线程协同完成，每个GPU线程计算结果矩阵的一个元素，然后加到变量sum上。4个GPU线程都计算完毕就会重复第三、四步，缓存
和
，计算
。如果还有其他子块，那么这个过程会一直重复下去。最终，每个GPU线程汇总矩阵
的一个元素。
代码 7
最后一步，把汇总的结果保存到矩阵C，如代码8所示。
代码 8
至此，我相信你已经深刻地体会到分组模式的复杂性。CPU拥有更强的控制部件和更大的缓存区域，可以预测和决定应该缓存哪些数据，而GPU则把原本属于它们的空间留给更多的运算部件，把缓存的控制权交给程序员，这意味着缓存的逻辑将会渗透到业务的逻辑，从而增加了代码的复杂性。
那么，这样做是否值得？我们可以算一下，在代码4里，av、bv和cv都是位于全局内存，每次访问都要1000个GPU时钟周期， 读取N次av和bv，写入一次cv，总共耗时2000 * N + 1000个GPU时钟周期，当N = 4时，总共耗时9000个GPU时钟周期。在代码6、7、8里，at和bt都是位于共享内存，每次访问只要10个GPU时钟周期，读取W次av和bv，写入W次at和bt，读取W * T次at和bt，写入一次cv，总共耗时 (2000 + 20 + 20 * T) * W + 1000，当N = 4，T = 2时，W = 2，总共耗时5120个GPU时钟周期，约为简单模式的56.89%，性能的改善非常明显。如果我们增加矩阵和子块的大小，这个差距就会更加明显，令N = 1024，T = 16，则简单模式总共耗时2049000个GPU时钟周期，而分组模式总共耗时150760个GPU时钟周期，后者是前者的7.36%。
当然，这并不是简单模式和分组模式在性能上的精确差距，因为我们还没考虑访问寄存器和算术运算的耗时，但这些操作的耗时和访问全局内存的相比简直就是小巫见大巫，即使把它们考虑进去也不会对结果造成太大影响。
你可能会问的问题
1. 什么时候不能使用分组模式？
分组模式最高只能处理三维的数据结构，四维或者更高的数据结构必须使用简单模式，事实上，简单模式会把四维或者更高的数据结构换算成三维的。如果算法没有反复从全局内存访问相同的数据，那也不必使用分组模式。
2. 分组的限制有哪些？
分组的大小不能超过1024，对应的维度必须能被整除，第一、二维度的大小不能超过1024，第三维度的大小不能超过64，分组的总数不能超过65535。
3. 分组的大小越大越好吗？
不是的，我们使用分组模式主要是为了使用共享内存，一般情况下，分组的大小和它使用的共享内存成正比，每个流多处理器的共享内存是有限的，比如说，NVIDIA的GK110 GPU的最大规格是48K，每个流多处理器最多可以同时容纳16个分组，这意味着每个分组最多只能使用3K，如果每个分组使用4K，那么每个流多处理器最多只能同时容纳12个分组，这意味着流多处理器的计算能力没被最大限度使用。
4. 在设定分组的大小时需要考虑Warp Size吗？
在计算域允许的情况下尽可能考虑，NVIDIA的Warp Size是32，AMD的Wavefront Size是64，因此，分组的大小最好是64的倍数。
*声明：本文已经首发于InfoQ中文站，版权所有，《遇见C++ AMP：GPU的线程模型和内存模型》，如需转载，请务必附带本声明，谢谢。