在本章中，我们将从一个相对简单的应用的背景和问题阐述开始，这个应用传统上一直受限于主流计算系统的有限能力。我们展示了并行执行不仅能够加速现有方法，还能让应用专家们采纳一种已知能带来益处，但由于计算需求过高而被此前忽略的方法。这种方法代表了一类日益重要的计算方法，它们能够从海量观测数据中推导出未知值的统计最优估计。我们将使用一个源自此类方法的示例算法及其实现源代码，来阐述开发者如何系统地确定内核并行结构、将变量分配到不同类型的内存中、规避硬件限制、验证结果，并评估性能提升所带来的影响。
17.1 背景
磁共振成像（MRI）是一种常用的医疗程序，用于安全、无创地探测人体所有区域的生物组织结构和功能。通过MRI生成的图像对临床和研究环境都产生了深远的影响。MRI包括两个阶段：采集（扫描）和重建。在采集阶段，扫描仪沿预定轨迹在k空间域（即空间频率域或傅里叶变换域）对数据进行采样。然后，在重建阶段，这些样本被转换为所需的图像。直观地说，重建阶段是根据扫描仪收集到的k空间观测数据来估计组织形状和纹理的过程。
MRI的应用常受限于高噪声水平、显著的成像伪影和/或较长的数据采集时间。在临床环境中，较短的扫描时间不仅可以提高扫描仪的吞吐量，还可以减轻患者的不适，并有助于减轻运动相关的伪影。高图像分辨率和保真度非常重要，因为它们能够实现病变的早期检测，从而改善患者的预后。然而，短扫描时间、高分辨率和高信噪比（SNR）的目标往往相互冲突；提升其中一项指标往往是以牺牲另外一项或两项指标为代价的。我们需要新的技术突破来实现所有这三个维度的同步改进。本研究提出了一个大规模并行计算提供了这种突破的案例。
MRI的物理原理请读者参考《Liang and Lauterbur (1999)》等MRI教科书。在本案例研究中，我们将重点关注重建阶段的计算复杂性，以及复杂性如何受到k空间采样轨迹的影响。MRI扫描仪使用的k空间采样轨迹可以显著影响重建图像的质量、重建算法的时间复杂性以及扫描仪采集样本所需的时间。式（17.1）显示了将k空间样本与重建图像关联起来的公式，适用于一类重建方法：
\[m(r) = \sum_j [ W(k_j) * s(k_j) * e^{i2 \pi k_j \cdot r} ] \;\;\; (17.1)\]在式（17.1）中，m(r)是重建图像，s(k)是测量的k空间数据，W(k)是权重函数，用于处理非均匀采样；即，W(k)会降低采样点密度较高的k空间区域数据的影响。对于这类重建，W(k)还可以作为一种窗函数（apodization filtering function），用于减小噪声的影响并减少由于有限采样导致的伪影。
如果在理想条件下，数据是在k空间中均匀间隔的笛卡尔网格点上采集的，那么W(k)权重函数是一个常数，因此可以从式（17.1）的求和中提出来。此外，对于均匀间隔的笛卡尔网格样本，式（17.1）中的指数项在k空间中是均匀间隔的。因此，m(r)的重建成为s(k)上的逆快速傅里叶变换（FFT），这是一种计算效率极高的方法。在这些均匀间隔的笛卡尔网格点上测量的数据集被称为笛卡尔扫描轨迹。图17.1A描绘了一个笛卡尔扫描轨迹。在实践中，笛卡尔扫描轨迹易于在扫描仪上实现，并被广泛应用于当今的临床环境中。
图17.1 扫描仪 k 空间轨迹及其相关的重建策略：(A) 笛卡尔轨迹配合 FFT 重建，(B) 螺旋（或一般非笛卡尔）轨迹后接网格化，以实现 FFT 重建，(C) 螺旋（非笛卡尔）轨迹配合基于线性求解器的重建。
尽管笛卡尔扫描数据的逆FFT重建在计算上非常高效，但非笛卡尔扫描轨迹通常具有优势，例如对患者运动的敏感性降低、更好地提供自校准场不均匀信息，以及对扫描仪硬件性能要求降低。因此，已经提出了诸如螺旋线（如图17.1C所示）、径向线（也称为投影成像）和玫瑰花形等非笛卡尔扫描轨迹，以减少运动相关的伪影并解决扫描仪硬件性能限制。这些改进最近使得重建图像的像素值能够用于测量诸如组织化学异常等细微现象，而这些现象在成为解剖病理学特征之前是难以发现的。
图17.2 非笛卡尔 k 空间样本轨迹和准确的基于线性求解器的重建，为令人兴奋的医疗应用带来了新的能力。
图17.2显示了这样一个基于MRI重建的测量，它生成了钠的图谱，钠是一种在正常人体组织中受到严格调控的物质。该信息可用于跟踪中风和癌症治疗过程中组织的健康状况。由于钠在人体组织中的含量远低于水分子，可靠地测量钠水平需要通过更多的样本来提高信噪比，因此需要利用非笛卡尔扫描轨迹来减轻额外的扫描时间。改进后的信噪比使得可靠地收集人体组织中化学物质（例如钠）的活体（in vivo）浓度数据成为可能。钠浓度的变化或移动表明疾病发展或组织坏死的早期迹象。例如，图17.2所示的人脑钠图可用于早期指示脑肿瘤组织对化疗方案的反应性，从而实现个体化医疗。
从非笛卡尔轨迹数据重建图像带来了挑战和机遇。主要的挑战在于指数项不再均匀间隔；求和形式不再是FFT。因此，不能再通过直接对k空间样本应用逆FFT来进行重建。在一种常用的方法网格化（gridding）中，样本首先被插值到均匀的笛卡尔网格上，然后使用FFT进行重建（参见图17.1B）。例如，一种卷积方法用于网格化，它获取一个k空间数据点，将其与一个网格化卷积掩模进行卷积，并将结果累加到笛卡尔网格上。正如我们在第7章所见，卷积的计算量相当大，是进行大规模并行计算的重要模式。读者已经具备了使用并行计算加速卷积网格化计算的技能，从而促进了当前FFT方法应用于非笛卡尔轨迹数据。
在本章中，我们将介绍一种迭代式、统计最优的图像重建方法，该方法可以准确地模拟成像物理并限定所得图像像素值中的噪声误差。随着大数据分析的兴起，此类统计最优方法正变得越来越重要。然而，与网格化相比，由于其过高的计算需求，此类迭代重建方法对于大规模三维（3D）问题一直不切实际。最近，由于GPU的广泛应用，这些重建方法在临床环境中变得可行。过去使用高端串行CPU需要花费数小时才能重建中等分辨率图像的迭代重建算法，现在使用CPU和GPU仅需数分钟，这一延迟在临床环境中是可以接受的。
17.2 迭代重建
Haldar 和 Liang 提出了一种基于线性求解器的迭代重建算法（Stone 等人，2008），用于非笛卡尔扫描数据，如图 17.1C 所示。该算法允许明确地对扫描仪数据采集过程的物理特性进行建模，从而可以减少重建图像中的伪影。然而，它的计算成本很高。我们以此为例，来说明那些创新但因计算时间过长而未被认为实用的方法。我们将展示大规模并行执行可以将重建时间缩短到分钟级别，从而使新的成像能力（例如钠成像）能够在临床环境中部署。
图17.3 一种用于重建非笛卡尔 k 空间样本数据的基于线性求解器的迭代方法。
图 17.3 展示了基于迭代线性求解器的重建方法的准贝叶斯（quasi-Bayesian）估计问题公式的解，其中 $\rho$ 是一个向量，包含重建图像的体素值；$F$ 是一个对成像过程的物理特性进行建模的矩阵；$D$ 是来自扫描仪的数据样本向量；$W$ 是一个可以包含解剖约束等先验信息的矩阵。$F^H$ 和 $W^H$ 分别是 $F$ 和 $W$ 的厄米特转置（或共轭转置），通过对矩阵进行转置后取每个元素的复共轭（$a + ib$ 的复共轭是 $a - ib$）得到。在临床环境中，$W$ 中表示的解剖约束来源于患者的一次或多次高分辨率、高信噪比的水分子扫描。这些水分子扫描揭示了诸如解剖结构位置等特征。矩阵 $W$ 就是从这些参考图像中推导出来的。问题是：给定所有其他矩阵求解向量 $\rho$。
从表面上看，图 17.3 中问题公式的计算解应该非常简单。它涉及矩阵乘法和加法（$F^H F + \lambda W^H W$）、矩阵-向量乘法（$F^H D$）、矩阵求逆 $((F^H F + \lambda W^H W)^{-1})$，最后是矩阵乘法 $((F^H F + \lambda W^H W)^{-1} F^H D)$。然而，这些矩阵的尺寸使得这种直截了当的方法非常耗时。$F^H$ 和 $F$ 矩阵的维度由 3D 重建图像中的体素数量和重建中使用的 k 空间样本数量决定。即使在一个适度的 $128^3$ 体素重建中，F 矩阵就有 $128^3 \approx 200$ 万列，每列有 $N$ 个元素，其中 $N$ 是使用的 k 空间样本数量（即 $D$ 的大小）。显然，$F$ 是一个极其巨大的矩阵。当试图使用迭代求解器方法来估计海量噪声观测数据的主要贡献因素时，通常会遇到如此巨大的维度，这在大数据分析中很常见。
由于涉及的矩阵尺寸如此之大，以至于使用高斯消元法等直接解法对图 17.3 中方程进行矩阵运算在实践中是无法处理的。因此，更倾向于采用迭代的矩阵求逆方法，例如共轭梯度（CG）算法。CG 算法通过迭代求解图 17.3 中的方程来重建图像中的 $\rho$。在每次迭代中，CG 算法更新当前的图像估计值 $\rho$，以改进准贝叶斯成本函数的值。CG 技术的计算效率主要由涉及 $F^H F + \lambda W^H W$ 和 $\rho$ 的矩阵-向量乘法运算的效率决定，因为这些运算在 CG 算法的每次迭代中都是必需的。
幸运的是，矩阵 $W$ 通常具有稀疏结构，这允许 $W^H W$ 的高效实现，而矩阵 $F^H F$ 是Toeplitz矩阵，可以通过 FFT 实现高效的矩阵-向量乘法。Stone 等人（2008）提出了一种 GPU 加速方法来计算 $Q$，这是一种数据结构，可以让我们在不实际计算 $F^H F$ 本身的情况下快速计算涉及 $F^H F$ 的矩阵-向量乘法。$Q$ 的计算在高端 CPU 核心上可能需要数天。由于 $F$ 对成像过程的物理特性进行建模，它只需要针对给定的扫描仪和规划的轨迹计算一次。因此，$Q$ 只需计算一次，可用于使用相同扫描轨迹的多次扫描。
计算 $F^H D$ 的矩阵-向量乘法所需时间大约比计算 $Q$ 少一个数量级，但对于一个 $128^3$ 体素的重建，在高端顺序 CPU 上仍可能需要大约 3 小时。回想一下，$D$ 是来自扫描仪的数据样本向量。因此，由于每次图像采集都需要计算 $F^H D$，理想情况下需要将 $F^H D$ 的计算时间缩短到几分钟。我们将展示这一过程的细节。事实证明，$Q$ 的核心计算结构与 $F^H D$ 的核心计算结构完全相同；只是 $Q$ 涉及矩阵乘法而非仅仅矩阵-向量乘法，因此计算量大得多。因此，从并行化的角度来看，我们只需讨论其中一个。我们将重点关注 $F^H D$，因为这是每次数据采集都需要运行的部分。
图 17.3 中的“find $\rho$”步骤执行基于 $F^H D$ 的实际 CG 算法。正如我们前面解释的，预先计算 $Q$ 使得这一步骤的计算强度比 $F^H D$ 小得多，在串行 CPU 上，它占每次图像重建执行时间的不到 1%。因此，在本章中，我们将把 CG 求解器排除在并行化范围之外，重点关注 $F^H D$。但是，我们将在本章末尾重新审视它的状态。
17.3 计算$F^H D$
图17.4 计算$F^H D$
图 17.4 展示了计算 $F^H D$ 核心步骤所需数据结构的顺序 C 语言实现。计算从一个外层循环开始，该循环遍历 k 空间样本（第 01 行）。
快速浏览图 17.4 可以发现，$F^H D$ 的 C 语言实现是加速的绝佳候选，因为它展现出大量数据并行性。
该算法首先在 k 空间的当前样本点计算 Mu 的实部和虚部（rMu 和 iMu）。然后进入一个内层 n 循环，计算当前 k 空间样本对图像空间中每个体素的 $F^H D$ 实部和虚部的贡献。请记住，$M$ 是 k 空间样本的总数，$N$ 是重建图像中的体素总数。$F^H D$ 在任何一个体素的值都取决于所有 k 空间样本点的值。然而，$F^H D$ 的任何体素元素都不依赖于 $F^H D$ 的任何其他体素元素。因此，$F^H D$ 的所有元素都可以并行计算。具体来说，外层循环的所有迭代都可以并行进行，内层循环的所有迭代也可以并行进行。但是，内层循环的计算依赖于外层循环同一迭代中先前语句完成的计算。
尽管该算法具有丰富的固有并行性，但潜在的性能瓶颈依然明显。
首先，在计算 $F^H D$ 元素的循环中，浮点运算与内存访问的比率最好情况下仅为 $0.75 \text{ OP/B}$，最差情况下为 $0.25 \text{ OP/B}$。最好情况假设 $\sin$ 和 $\cos$ 三角函数运算是使用分别需要 13 个和 12 个浮点运算的五项泰勒级数计算的。最差情况假设每个三角函数运算在硬件中作为单个操作计算。正如我们在第 5 章《内存架构和数据局部性》中所见，要使内核不受到内存带宽的限制，需要一个高得多的浮点算术与全局内存访问比率。因此，除非该比率得到大幅提高，否则内存访问将明显限制内核的性能。
其次，浮点算术与浮点三角函数之比仅为 $13:2$。因此，基于 GPU 的实现必须容忍或避免由于 $\sin$ 和 $\cos$ 运算的长延迟和低吞吐量而导致的停顿。如果没有好的方法来降低三角函数的成本，性能很可能会被花费在这些函数上的时间所主导。
我们现在已经准备好将 $F^HD$ 从顺序 C 代码转换为 CUDA 内核 的各个步骤。
第一步：确定内核的并行结构
将图 17.4 中的循环转换为 CUDA 内核在概念上是直接明了的。由于图 17.4 的外层循环的所有迭代都可以并行执行，我们可以简单地将外层循环映射到 CUDA 线程上，将其转换为一个 CUDA 内核。图 17.5 展示了这种直接转换后的内核。
图17.5 $F^H D$内核的第一个版本
在该实现中，每个线程对应原外层循环的一次迭代，也就是说，每个线程负责计算一个 k 空间采样点（k-space sample） 对所有 F H D 元素的贡献。原外层循环共有 M 次迭代，而 M 的值可能达到数百万。因此，我们显然需要大量线程块（thread blocks），以产生足够多的线程来执行所有这些迭代。
为了方便性能调优，我们定义了一个常量 FHD_THREADS_PER_BLOCK，用于指定在调用 cmpFhD 内核时每个线程块中线程的数量。因此，在调用内核时，我们将使用：
- 网格大小（grid size） = M / FHD_THREADS_PER_BLOCK
- 块大小（block size） = FHD_THREADS_PER_BLOCK
在内核内部，每个线程通过以下常见公式计算它所负责的原外层循环迭代索引：
blockIdx.x * FHD_THREADS_PER_BLOCK + threadIdx.x
例如，假设有 1,000,000 个 k 空间采样点，我们决定每个块使用 1024 个线程，则网格大小为 1,000,000 / 1024 = 977 个块，块大小为 1024。每个线程的 m 值计算公式为：
m = blockIdx.x * 1024 + threadIdx.x
虽然图 17.5 的内核能充分利用并行性，但它存在一个主要问题：所有线程都会向所有 rFhD 和 iFhD 体素（voxel）元素写入数据。
这意味着，为了防止线程之间在更新体素值时互相覆盖（即“相互踩踏”），必须在内循环中使用全局内存上的原子操作（atomic operations）（见第 10–11 行）。
正如我们在第 9 章《并行直方图》中看到的那样，在全局内存上大量使用原子操作会严重降低并行执行的性能。
此外，rFhD 和 iFhD 数组的大小（即重建图像的体素总数）也使得在共享内存中进行私有化（privatization）不可行，因此我们需要探索其他方案。
在这种计算模式下，每个线程从一个输入元素（在本例中是一个 k 空间采样点）出发，更新所有或多个输出元素（重建图像中的体素），这种方式称为散射方法（scatter approach）。
直观地说，每个线程将一个输入的影响“散射”到多个输出值上。 然而，这种方法的缺点是：多个线程可能同时更新相同的输出元素，从而导致相互覆盖，因此必须使用原子操作来保证正确性，而这会显著影响并行性能。
一种更优的替代方案是 聚合方法（gather approach）。 在这种方式中，每个线程计算一个输出元素，通过汇总所有输入元素的贡献来获得最终值。 这样，每个线程只更新自己对应的输出元素，不会与其他线程发生冲突。
在我们的应用中，基于聚合方法的并行化意味着：
- 每个线程负责计算一个
rFhD/iFhD元素对； - 每个线程从所有 k 空间采样点收集贡献。
这样就避免了线程之间的干扰，也不再需要原子操作。
要采用聚合方法，我们需要使 n 循环 成为外层循环，以便将其每次迭代分配给一个线程。
这可以通过交换内外层循环（loop interchange）来实现。
交换后，每个新的外层循环迭代处理一个 rFhD / iFhD 元素对，而新的内层循环则累积所有 k 空间采样点对该元素对的贡献。
这种结构转换称为 循环交换（loop interchange）。
它要求循环是“完美嵌套的”（perfectly nested），即外层 for 语句和内层 for 语句之间不能有其他语句。
然而，在图 17.4 的 F H D 代码中并非如此，因此我们需要先将 rMu 和 iMu 元素的计算移出循环结构。
从对图 17.4 的快速检查可以看出，可以使用一种称为循环分裂（loop fission）或循环拆分（loop splitting）的技术将 F H D 计算拆分为两个独立的循环，如图 17.6 所示。
图17.6 循环裂变 (Loop Fission) 应用于$F^H D$计算。
循环分裂的基本思想是：将一个循环的循环体分成两部分。在 F H D 的情况下，外层循环包含两个部分：
- 内层循环之前的语句；
- 内层循环本身。
通过循环分裂，我们可以将前者放入一个新的循环（第一阶段），后者放入另一个新的循环（第二阶段）。
需要注意的是，循环分裂会改变原循环中语句的执行顺序。在原循环中，每次迭代的两部分顺序执行；而分裂后，所有迭代的第一部分会先执行完，然后再执行所有迭代的第二部分。 读者可以验证，这种顺序变化并不会影响 F H D 的执行结果，因为每次迭代的第一部分并不依赖于任何先前迭代的第二部分结果。
循环分裂是高级编译器常用的优化转换，它依赖于对循环间依赖关系的自动分析。
经过循环分裂后，$F^HD$ 的计算可分为两个阶段：
- 第一阶段：单层循环，用于计算
rMu和iMu元素（供第二阶段使用）； - 第二阶段：计算
F H D元素，使用第一阶段计算的rMu和iMu。
每个阶段都可以分别转换为一个 CUDA 内核，并顺序执行（第二个内核依赖第一个的结果）。 由于两者在逻辑上本就需要顺序执行，因此这种拆分不会损失任何并行性。
图17.7 cmpMu内核
图 17.7 所示的 cmpMu() 内核实现了第一阶段的循环。
从顺序 C 代码到 CUDA 内核的转换相对简单：每个线程执行原始循环的一次迭代。
由于 M（k 空间采样点数）可能非常大，我们需要使用多个线程块。
假设每个线程块最多可包含 1024 个线程，则可以设置常量 MU_THREADS_PER_BLOCK = 1024，并使用：
- 块大小（block size） =
MU_THREADS_PER_BLOCK - 网格大小（grid size） =
M / MU_THREADS_PER_BLOCK
例如，当有 1,000,000 个采样点时，可配置为 1024 线程每块、977 个块（1,000,000 / 1024 ≈ 977）。
在内核中，每个线程可根据其 blockIdx 和 threadIdx 计算自己对应的原始循环迭代号：
iteration = blockIdx.x * MU_THREADS_PER_BLOCK + threadIdx.x
例如，当 MU_THREADS_PER_BLOCK = 1024 时：
- 线程
(blockIdx.x=0, threadIdx.x=37)对应第 37 次迭代； - 线程
(blockIdx.x=5, threadIdx.x=2)对应第 5122 次迭代（5×1024+2）。
这样，通过使用该迭代号访问 Mu、Phi 和 D 数组，可以确保这些数组的访问模式与原 C 代码循环一致。
由于每个线程只写入自己的 Mu 元素，不会发生写冲突。
确定第二个核函数的结构需要多做一些工作。查看图 17.6 中的第二个循环可以发现，设计第二个核函数至少有三种选择。 在第一种选择中，每个线程对应内层循环的一次迭代。这种方案创建的线程数量最多，因此可以利用最大程度的并行性。然而，线程总数将达到 N×M，其中 N 在数百万级，M 在几十万级。它们的乘积会导致网格中出现过多的线程，远超出充分利用设备所需的数量。
第二种选择是让每个线程实现外层循环的一次迭代。与第一种方案相比，这种方案使用的线程更少。该方案生成 M 个线程，而不是 N×M 个线程。由于 M 对应于 k 空间采样点的数量，而通常用于计算 FHD 的采样点数很大，大约在十万量级，因此这种方案仍然能发挥较高程度的并行性。然而，这个核函数会遇到与图 17.5 中核函数相同的问题。也就是说，每个线程都会写入所有的 rFhD 和 iFhD 元素，从而造成线程之间极多的写入冲突。正如图 17.5 中的情况，图 17.8 的代码需要使用原子操作，这将显著降低并行执行的速度。因此，这种方案的效果并不好。
图17.8 $F^HD$内核的第二种方法
第三种选择是让每个线程计算一对 rFhD 和 iFhD 输出元素。这种方案要求我们交换内外层循环的位置，然后让每个线程实现新的外层循环的一次迭代。该变换如图 17.9 所示。循环交换是必要的，因为由 CUDA 线程实现的循环必须是外层循环。经过循环交换后，每次新的外层循环迭代都将处理一对 rFhD 和 iFhD 输出元素，而每次内层循环则会累加所有输入元素对这对输出元素的贡献。
图17.9 $F^HD$计算的循环交换
在这里进行循环交换是允许的，因为两个层级的所有循环迭代之间都是相互独立的。它们可以以任意顺序执行。循环交换仅改变迭代的执行顺序，只要这些迭代之间无依赖关系，这样的操作就是合法的。通过这种选择，我们可以将新的外层循环转换为一个由 N 个线程执行的核函数。由于 N 对应重建图像中的体素数，对于高分辨率图像而言，N 的值可能非常大。例如，对于一个 128³ 的图像，共有 128³ = 2,097,152 个线程，可实现高度的并行性。对于更高的分辨率（如 512³），我们可能需要使用多维网格、启动多个网格，或让单个线程处理多个体素。
在第三种方案中，每个线程都只对自己的 rFhD 和 iFhD 元素进行累加，因为每个线程都有唯一的 n 值，因此线程之间不会发生冲突。这使得第三种方案成为三种方案中最优的选择。
图17.10 $F^HD$内核的第三种方法
由交换循环得到的核函数如图 17.10 所示。外层循环已被去除；每个线程负责外层（n）循环的一次迭代，其中 n = blockIdx.x × FHD_THREADS_PER_BLOCK + threadIdx.x。一旦确定该迭代的 n 值，线程便根据该 n 值执行内层（m）循环。
该核函数可通过指定全局常量 FHD_THREADS_PER_BLOCK 来设置每个块的线程数量。假设变量 N 存储重建图像的体素数，那么 N / FHD_THREADS_PER_BLOCK 个块就能覆盖原始循环的所有 N 次迭代。
例如，如果有 2,097,152 个体素，可以使用每块 1024 个线程、共 2,097,152 / 1024 = 2048 个块的配置来调用核函数。在图 17.10 中，这通过将 FHD_THREADS_PER_BLOCK 设为 1024，并在核函数调用时将其作为块大小、将 N / FHD_THREADS_PER_BLOCK 作为网格大小来实现。
第二步：突破内存带宽限制
图 17.10 所示的简单 cmpFhD 核函数的性能将明显优于图 17.5 和图 17.8 中的核函数，但由于受到内存带宽的限制，其加速效果仍然有限。快速分析表明，其执行受限于每个线程较低的计算与全局内存访问比。
在原始循环中，每次迭代至少会进行 14 次内存访问：kx[m]、ky[m]、kz[m]、x[n]、y[n]、z[n]、rMu[m]（两次）、iMu[m]（两次）、rFhD[n]（读写各一次）以及 iFhD[n]（读写各一次）。与此同时，每次迭代执行大约 13 次浮点乘法、加法或三角函数运算。因此计算与全局内存访问比为 13 / (14×4) = 0.23 OP/B，这个值太低（根据第 5 章“内存结构与数据局部性”的分析）。
我们可以通过将部分数组元素分配给自动变量来立即提高计算与全局内存访问比。正如第 5 章所述，自动变量将存放在寄存器中，从而将全局内存的读写转化为片上寄存器的读写。快速回顾图 17.10 的核函数可知，对于每个线程，x[n]、y[n] 和 z[n] 元素在循环的所有迭代中都保持不变（第 05–06 行）。这意味着我们可以在进入循环前将这些元素加载到自动变量中。然后，核函数在循环中使用这些寄存器变量，从而将全局内存访问转化为寄存器访问。
此外，循环中会反复读取和写入 rFhD[n] 和 iFhD[n]。我们可以让循环的迭代仅访问两个自动变量，并在循环结束后再将这些变量的内容写回到 rFhD[n] 和 iFhD[n]。这样得到的代码如图 17.11 所示。通过为每个线程多使用 5 个寄存器，我们将每次迭代的内存访问从 14 次减少到 7 次，使计算与全局内存访问比从 0.23 OP/B 提高到 0.46 OP/B。这是一个不错的改进，也体现了寄存器资源的合理利用。
图17.11 使用寄存器减少$F^HD$内核的全局内存访问
需要注意的是，寄存器使用量可能会限制占用率（即每个流式多处理器 SM 上能同时运行的线程块数量）。在核函数中每个线程增加 5 个寄存器使用量，相当于每个线程块增加 5×FHD_THREADS_PER_BLOCK 个寄存器。假设每个块有 1024 个线程，则每个块增加了 5120 个寄存器使用量。由于每个 SM（版本 3.5 或更高）可容纳所有线程块共计 65,536 个寄存器，因此我们必须小心，进一步增加寄存器使用量可能会限制分配给 SM 的块数量。幸运的是，对于该核函数而言，寄存器使用量并不会成为并行性的限制因素。
我们希望进一步通过减少更多全局内存访问来提高 cmpFhD 核函数的计算与内存访问比。接下来的候选对象是 k 空间采样点 kx[m]、ky[m] 和 kz[m]。这些数组元素的访问方式与 x[n]、y[n] 和 z[n] 不同：在图 17.11 的循环中，每次迭代访问的 kx、ky 和 kz 元素都不同。这意味着我们无法将 k 空间元素加载到寄存器后在整个循环中反复使用，因此寄存器在此无效。
然而，我们应注意到 k 空间元素不会在核函数中被修改，而且每个元素会被网格中所有线程使用。这表明我们可以将这些 k 空间元素放入常量内存中，利用常量缓存来消除大多数 DRAM 访问。
对图 17.11 的循环进行分析可知，k 空间元素确实非常适合存放在常量内存中。用于访问 kx、ky 和 kz 的索引是 m，而 m 与 threadIdx 无关，这意味着同一个 warp 中的所有线程都会访问相同的 kx、ky、kz 元素。这正是常量缓存理想的访问模式：每当一个元素被加载到缓存中，它至少会被当前 warp 的 32 个线程全部使用。因此，每 32 次常量内存访问中，至少有 31 次可由缓存满足，相当于消除了 96% 以上的全局内存访问。更好的是，从缓存中读取的常量可以广播给 warp 中所有线程，使得访问常量内存的效率几乎与访问寄存器相当。
不过，将 k 空间元素放入常量内存还存在技术问题：常量内存容量仅为 64KB，而 k 空间样本的大小可能大得多，往往达到数百万级。常见的解决办法是将大数据集分块，每块大小不超过 64KB。开发者需要重新组织核函数，使其被多次调用，每次只处理一块数据。这对 cmpFhD 核函数来说很容易实现。
仔细观察图 17.11 的循环可以发现，所有线程都会顺序遍历 k 空间采样数组，即在每次迭代中，网格中的所有线程都访问同一个 k 空间元素。对于大型数据集，循环只需执行更多次。因此，我们可以将循环划分为多个部分，每部分处理一块可放入 64KB 常量内存的数据。主机代码将多次调用核函数，每次调用前通过 cudaMemcpyToSymbol() 将新的数据块传入常量内存，如图 17.12 所示。（在较新的设备和 CUDA 版本中，可以通过 const __restrict__ 声明内核参数，使其数据自动驻留于只读缓存，从而达到与常量内存类似的效果。）
图17.12 将 k 空间数据分块以适配常量内存的主机端代码序列。
在图 17.12 中，cmpFhD 核函数在一个循环中被调用。代码假定 kx、ky、kz 数组存放在主机内存中，其维度为 M。每次迭代，主机代码调用 cudaMemcpyToSymbol() 将一块 k 空间数据传入设备常量内存（详见第 7 章“卷积”），然后调用核函数处理该块数据。当 M 不是 CHUNK_SIZE 的整数倍时，主机代码需要额外进行一次数据传输和一次核函数调用来处理剩余数据。
图17.13 使用常量内存的修订版$F^HD$内核
图 17.13 展示了一个从常量内存访问 k 空间数据的改进核函数。注意此时核函数参数列表中不再包含 kx、ky、kz 指针。kx_c、ky_c、kz_c 数组作为全局变量用 __constant__ 关键字声明（见图 17.12）。通过从常量缓存访问这些元素，核函数现在仅需对 rMu 和 iMu 数组进行四次全局内存访问。编译器通常会识别出这些访问实际上只涉及两个位置，即各对 rMu[m] 和 iMu[m] 各一次访问，其值被存入临时寄存器变量供后续使用。最终内存访问次数减少为 2 次，计算与内存访问比提升至 1.63 OP/B。这虽然仍不算理想，但足以使内存带宽不再成为性能的唯一瓶颈。稍后我们还会看到更多可提高效率的优化。
如果我们在某些设备上运行图 17.12 与 17.13 的代码，会发现性能提升不如预期。这是因为常量缓存的性能未达到预期，问题源于常量缓存的设计与 k 空间数据的内存布局。
图17.14 k 空间数据布局对常量缓存效率的影响：（A）将 k 空间数据存储在独立数组中；（B）将 k 空间数据存储在由结构体组成的数组中。
如图 17.14A 所示，每个常量缓存行用于存储多个连续字。这样的设计可降低硬件成本。当某个元素被加载到缓存中时，其周围的若干元素也会一并被加载（图中阴影部分表示）。在每个 warp 的一次迭代中，通常需要三个缓存行来保证高效执行。
在典型的执行过程中，一个 SM 上会同时运行大量 warp。由于不同 warp 所处的迭代阶段不同，可能需要同时占用许多缓存行。例如，若每个线程块包含 1024 个线程，并在每个 SM 上同时运行两个块，则共有 (1024 / 32) × 2 = 64 个 warp 并发执行。如果每个 warp 至少需要 3 个缓存行维持执行效率，则最坏情况下需要 64 × 3 = 192 个缓存行。即使假设平均每 3 个 warp 处于相同迭代阶段、可共享缓存行，也仍需约 64 个缓存行——这就是所谓的活动 warp 工作集。
由于成本限制，一些设备的常量缓存仅包含 32 条缓存行。当缓存行不足以容纳整个工作集时，不同 warp 访问的数据会相互竞争缓存空间。等到某个 warp 进入下一次迭代时，所需数据可能已被其他 warp 的数据替换出去。结果，某些设备的常量缓存容量确实不足以同时容纳所有活动 warp 的数据，从而未能有效减少全局内存访问。
针对缓存使用效率低的问题，已有研究提出可通过调整 k 空间数据的内存布局来解决。解决方案如图 17.14B 所示，对应代码见图 17.15 和 17.16。其核心思想是：不再将 k 空间的 x、y、z 分量分别存放在三个数组中，而是将它们封装为结构体的成员，再组成一个数组。这种声明方式称为结构体数组（array of structures）。
图17.15 调整 k 空间数据布局以提升缓存效率。
图 17.15（第 01–03 行）展示了数组的声明形式。假设内存分配与初始化代码（未显示）已正确地将 k 空间数据的 x、y、z 分量填充到结构体字段中。通过将这三个分量存放在连续的常量内存位置中，每次 warp 迭代所需的三个分量可同时存入一个缓存行，从而减少支持所有活动 warp 执行所需的缓存行数量。
由于现在只需一个数组保存所有 k 空间数据，因此主机只需一次 cudaMemcpyToSymbol() 调用即可将整个数据块拷贝至常量内存。假设每个 k 空间样本为单精度浮点数，则传输大小从原来的 4×CHUNK_SIZE 调整为 12×CHUNK_SIZE，以反映三个分量同时传输。
采用新数据结构布局后，核函数也需相应修改，按新格式访问数据。新的核函数如图 17.16 所示，其中 kx[m] 改为 k[m].x，ky[m] 改为 k[m].y，以此类推。此处虽是细微修改，但可显著提升部分设备的执行速度。
图17.16 在$F^HD$内核中调整以适应 k 空间数据的内存布局。
第三步：使用硬件三角函数
CUDA 提供了数学函数的硬件实现，其吞吐量远高于对应的软件实现。GPU 提供这些三角函数（如 sin() 和 cos()）的硬件实现的动机之一，是为了加快图形应用中的视角变换计算。这些函数作为硬件指令由 SFU（特殊功能单元，Special Function Units）执行。使用这些函数的过程相当简单。
在 cmpFhD 内核的情况下，我们所需做的就是将对 sin() 和 cos() 函数的调用替换为它们的硬件版本：__sin() 和 __cos()（函数名前加两个下划线 "__"）。这些是编译器可识别的内建函数（intrinsic functions），会被翻译成 SFU 指令。
由于这些函数被调用于高频执行的循环体中，我们预期此更改将带来显著的性能提升。修改后的 cmpFhD 内核如图 17.17 所示。
图17.17 使用硬件函数 __sin() 和 __cos()。
然而，我们需要注意的是，从软件函数切换到硬件函数会降低精度。当前的硬件实现精度低于软件库（详细信息见《CUDA C 编程指南》）。对于 MRI（磁共振成像）应用，我们必须确保硬件实现的精度足够，如图 17.18 所定义。测试过程涉及一个虚拟物体的“理想”图像（I₀），有时称为幻影物体（phantom object）。
图17.18 用于验证硬件函数精度的度量指标。I₀ 表示理想图像，I 表示重建图像，PSNR 表示峰值信噪比。
我们使用逆向过程生成相应的“扫描”k 空间数据（合成数据）。该合成的扫描数据再由所提出的重建系统处理，以生成重建图像（I）。理想图像与重建图像中各体素的值被输入到图 17.18 的峰值信噪比（PSNR）公式中。
测试的通过标准取决于图像的具体应用。在我们的实验中，我们与临床 MRI 专家合作，确保因使用硬件函数而导致的 PSNR 变化仍在其应用所接受的范围之内。 在医生利用图像进行伤情判断或疾病评估的应用中，还需要对图像质量进行目视检查。图 17.19 展示了原始“真实”图像的视觉对比结果。可以看到，CPU 双精度与单精度实现的 PSNR 都为 27.6 dB，远高于该应用的可接受水平。目视检查同样表明，重建图像与原始图像高度一致。
图 17.19 还清晰显示了迭代重建相较于简单双线性插值 gridding/iFFT 的优势。使用简单 gridding/iFFT 重建的图像 PSNR 仅为 16.8 dB，明显低于迭代重建方法的 27.6 dB。图 17.19（图像 2）中可见的严重伪影，会显著影响图像在诊断中的可用性，而这些伪影在迭代重建图像中并不存在。
当我们从 CPU 的双精度算术切换到单精度算术时，PSNR 未出现明显下降，仍保持在 27.6 dB。当我们将三角函数从软件库切换到硬件单元时，PSNR 仅轻微下降，从 27.6 dB 降至 27.5 dB。这种细微的损失在应用可接受范围内。目视检查也确认了重建图像与原始图像相比未出现明显伪影。
步骤 4：实验性性能调优
到目前为止，我们尚未确定内核配置参数的最佳取值。其中一个配置参数是每个线程块的线程数。选择合适的线程数量对于充分利用每个 SM 的线程容量至关重要。另一个配置参数是for 循环体的展开（unroll）次数，如图 17.17 第 07 行所示。
我们可以通过指令 #pragma unroll 后跟展开次数，来让编译器在循环上进行展开。
一方面，循环展开可以减少指令开销，并可能降低处理每个 k 空间采样数据所需的时钟周期数；另一方面，过度展开可能增加寄存器使用量，从而减少可在每个 SM 中并行容纳的线程块数。
需要注意的是，这些配置参数之间的影响并非相互独立。增加一个参数值，可能会占用本可用于增加另一个参数的资源。因此，必须以实验的方式联合评估这些参数。潜在的组合数可能非常多。
在$F^HD$的情况下，通过系统地搜索所有参数组合并选择运行时间最优的配置，与仅依据经验性趋势进行启发式调优相比，性能提升约 20%。 Ryoo 等人（2008）提出了一种基于帕累托最优曲线（Pareto optimal curve）的方法，用于筛除大部分性能较差的组合。
17.4 小结
本章介绍了将一个循环密集型应用（MRI 图像的迭代重建）从串行形式并行化和优化的关键步骤。 我们首先探讨了并行化的组织方式——散射（scatter）方法与聚合（gather）方法。我们展示了将散射方法转换为聚合方法是避免原子操作的关键，而原子操作会显著降低并行执行性能。 随后，我们讨论了实现聚合方法所需的实际技术，如循环分裂（loop fission）与循环互换（loop interchange）。
接着，我们介绍了多种优化技术的应用，包括：
- 将数组元素提升到寄存器中；
- 使用常量内存/缓存存放输入数据；
- 利用硬件函数提高并行内核性能。
通过这些优化，程序的性能提升了约 10³ 倍。
在并行化和优化之前，$F^HD$占据了几乎 100% 的执行时间。一个有趣的现象是：在优化后，共轭梯度（CG）求解器（图 17.3 中的 “求 ρ” 步骤）所耗时间反而超过了 $F^HD$。 这表明我们已大幅加速了 $F^HD$，要想进一步提速，就必须加速 CG 求解器。
成功并行化与优化后，$F^HD$仅占约 50% 的执行时间，另一半主要耗费在 CG 求解器上。这是实际应用并行化中的常见现象：当一些耗时阶段被成功加速后，原本次要的阶段就会成为新的性能瓶颈。