矢量化——又一次变得重要

下载Intel® Parallel Studio的免费试用版

Robert H. Dodds Jr., 名誉教授, 伊利诺伊大学

向量化曾是高性能计算的支柱，如今随着每个处理器核心上的SIMD（向量）寄存器宽度不断增加，它重新焕发了新的重要性：AVX2支持四个双精度浮点数，AVX-512支持八个。向量长度的增加为提高某些类别的软件性能提供了机会（例如，大规模有限元分析，用于模拟机械系统的非线性行为）。开源研究代码WARP3D（也在GitHub上）专注于此类模拟，以支持能源生产及相关行业中关键部件的安全性和寿命预测方法学的开发。其学科重点（3D非线性实体）、开源特性和可控的代码大小（1400个例程）吸引了工业界和学术界的研究人员。

WARP3D反映了当前的并行计算层次结构：使用MPI*访问集群中的多个节点，使用OpenMP*在每个节点上使用线程，以及在每个线程内进行向量化。用于预测复杂行为（例如，大型部件焊接产生的内部应力）的现实模拟需要在当前计算机上花费数小时的并行计算。为获得解决方案而进行的普遍存在的矩阵运算非常适合这种并行计算层次结构。在核心级别增加向量化使用量，可以对并发执行的线程和MPI秩产生直接的性能提升效果。

本文简要介绍了WARP3D中实现的向量化。该代码在过去20年中不断发展，采用了Fortran中许多不断变化的特性。我们使用了Intel® Parallel Studio XE（包括Intel® Fortran Compiler、Intel® Math Kernel Library和Intel® MPI Library）的组件来构建可执行文件，并随Linux*、Windows*和MacOS*的开源发行版一起提供。WARP3D采用非线性有限元法的隐式公式，这需要求解非常大且不断变化的线性方程组——Intel Math Kernel Library中的PARDISO*和CPARDISO*软件包提供了卓越的性能。此处描述的示例也展示了新的Intel® Advisor 2017的线形图分析功能，该功能能够更方便、更深入地探索循环级别的性能。

图1展示了一个法兰的三维有限元模型分解为多个单元域（每个域/MPI秩），并将域内的单元分配给块。模型中的每个离散单元都关联着数量不等且通常较小的永久和临时数组（例如，6x6、24x24、60x60、6x60），在中到大型模型求解过程中会产生数百万个此类数组。

图2A显示了一种传统的、多层次的数组结构（AoS）用于存储单元数据，单个单元的各种数组很可能在内存中是连续的。虽然概念上很方便，但相对较小的单元级矩阵限制了内部循环的向量长度，并且在处理所有单元的相同数组类型时，无法调整以适应缓存大小的影响。

如图2B所示，另一种块状数据结构（针对单元刚度矩阵K_e和B_e作为示例）定义了一种数组结构（SoA）方案，使用密集二维和三维数组，其首个维度由分配给块的单元数量决定，在WARP3D中称为span。因此，一个块中所有单元的K_e矩阵（以及单独的B_e矩阵）在内存中是连续的。

通过一个主列表设置各域的块所有权，以下序列说明了WARP3D中MPI与线程并行结合的概念

c$OMP PARALLEL DO PRIVATE( blk )
  do blk = 1, number_ele_blks
     if( blk_to_domain_map(blk) .ne. myrank ) cycle
     call process_a_blk( blk )
  end do
c$OMP END PARALLEL DO

WARP3D采用简单的管理器-工作者设计，通过MPI在秩之间同步任务。在上文中，管理器（在秩0上）确保所有工作者执行OpenMP并行循环，以处理它们所拥有的块。process_a_blk例程根据需要动态分配块状数组的工作版本到一个派生类型中，该类型声明为该例程的局部变量，因此在每个线程中都是局部的。构建/更新单元矩阵的通常复杂的多层次计算例程使用向量式编码——这意味着它们不知道OpenMP和MPI运行时环境。此处隐含的单元基础计算的数据独立性代表了有限元方法的一个关键特性。全局数组被线程只读访问；少数一个向量上的块收集操作指定$OMP ATOMIC UPDATE进行同步。WARP3D中的所有线程操作都通过上面所示的十几个循环/指令完成。

Cray*向量机的发展促使了上述块状方法在最初由美国国家实验室创建的有限元代码中的开发。随着这些向量机的消亡，代码演变为使用线程并行处理单元块，然后跨MPI秩进行处理，同时保留块状作为调整各种缓存架构影响的一种手段。

此处的一个示例说明了传统AoS和块状SoA方法的关键特性和性能。一系列矩阵乘法更新模型中每个单元的数组。对于单个单元，在此例中，B_e、D_e和K_e矩阵的大小分别为6x60、6x6、60x60（均为双精度浮点数且为密集矩阵）。图2B显示了块状存储方案，其中大型数组是Fortran派生类型的可分配组件。在具有不同单元类型（因此矩阵大小不同）的模型中，WARP3D创建同质块，所有单元类型相同。对K_e的更新形式为：K_e ¬ K_e + transpose(B_e)D_eB_e。在法兰模拟过程中，此更新会发生3亿到10亿次（D_e、B_e在每次更新时也会改变）。在此示例中，D_e 和K_e是对称的，因此计算仅限于下三角项。

图3A显示了传统AoS方案的核心代码。计算外层循环顺序遍历模型中的所有单元，调用子例程的不同版本以处理对称、非对称等情况，而子例程中的内层循环遍历单个单元矩阵的尺寸要小得多的维度。Fortran内置的matmul函数在第19行执行第一个矩阵乘法，内层循环遍历相对较小的单元数组来计算K_e下三角项的最终更新。编译器在第23至25行展开k=1,6内层循环，并根据第19行发出的注释优化报告中关于循环交换的大量消息，似乎将matmul替换为等效的源代码进行编译。[编辑注：请参阅The Parallel Universe第27期中的“使用Intel® AVX-512提高性能的向量化机会”，了解生成和解释编译器优化报告的建议。]

图3B中块状SoA方案的代码按以下顺序运行计算循环

1. 遍历所有块（未显示）
2. 遍历单个单元矩阵的维度
3. 遍历块中单元数量（span）的最内层循环

第10至28行的第一组循环在提供的内存区域中计算乘积db = D_eB_e，适用于块中的所有单元。直接实现需要四级循环结构。此处，最内层的两个循环进行了手动展开，本质上强制编译器向量化通常更长的i=1, span循环。对这些循环探索的其他排列会增加运行时间；编译器不知道每个循环的迭代次数，这限制了循环重新排序的有效性。第30至42行的第二组循环执行更新K_e ¬ K_e + transpose(B_e)db。在此，utsz=1820，使用索引向量icp将单元的下三角项映射到块状K_e的列。同样，最内层循环已手动展开以进行向量化i=1,span。

此例程中的两组循环（第10至28行和第30至42行）都具有单位步长数组引用和长（span）长度，为向量化而设置，这是WARP3D代码中普遍存在的特征。这些通常是WARP3D中单元级处理最耗时的循环；多年来，随着编译器优化和处理器功能的演进，循环和数组组织的多次重新设计已经发生。

图4总结了两种方法在单线程上运行时相对执行性能，AoS方案赋值为单位值。括号中的数字表示计算例程中数组的大小。相关的编译器选项为：ifort (17.0.2) 配合 -O3 –xHOST –align array64byte –ftz。计算机配备双Intel® Xeon® E5-2698 v3处理器（带AVX2），测试期间没有其他用户软件运行。该处理器具有缓存大小 L1 = 32 KB/核，L2 = 256 KB/核，L3 = 40 MB共享。驱动代码和计算例程位于不同的文件中，以防止内联效果。运行时间包括调用计算例程的努力，但不包括驱动程序设置时间。驱动程序为AoS和SoA测试分配了数千个数组集。每次调用时，驱动程序将一个随机选择的集传递给计算例程，所有情况下处理3200万个单元。多次运行，每次耗时2至4分钟，显示测量的CPU时间变异性为1%至2%。

AoS方案只需要一个结果；没有可调参数。SoA方案中每块的单元数（span）会变化，并伴有补偿的块数，以维持每种情况下处理的单元数量相同。SoA结果显示了块大小对向量处理长度和L1、L2、L3缓存利用率相互作用的预期强烈影响。当块大小为64和32个单元时，运行时间达到AoS时间的37%。当块大小为16个单元时，运行时间开始增加。对于所有块大小，数组都适合L3缓存；AoS中尺寸小得多的单单元数组几乎适合L2数据缓存。

图5显示了Intel® Advisor 2017执行的线形图分析的输出，并将关键情况的结果手动叠加到单个图上进行比较。注意两个轴都是对数刻度。Intel® Advisor中实现的带有图形工具的线形图分析提供了一种方便的机制，用于筛选大型代码中关键循环的性能（注意右上角的汉堡菜单）。将指针移到符号上会弹出显示源代码中循环位置、GFLOPS和每L1缓存字节的浮点运算次数的弹出窗口。单击符号会调出带有额外性能信息（例如，是否执行了循环的剥离循环和剩余循环）的循环源代码。

在K_e更新的内层循环（第33至41行）中，对于块大小（span）= 64，最佳性能为8.08 GFLOPS，该性能有效地处于L2缓存的线形图上，但远低于作为该循环核心的熔合乘加（FMA）向量指令的峰值速率。D_eB_e计算的内层循环（第11至27行）在块大小从1024减小到128时获得了惊人的3.7倍性能提升，但对于更小的块没有进一步提升。运行时间最长的循环（第33至41行）在块大小从1024减小到64时获得了2倍的性能提升。

正如计时结果（图4）所示，AoS的GFLOPS和缓存利用率（由蓝色填充和轮廓方形表示）小于块状SoA。用直接的ek = ek + matmul(transpose(b),matmul(d,b))替换例程中的代码，导致运行时间比图3A中的版本增加了30%。

在线形图上显示的8.08 GFLOPS的最佳性能相对于双精度FMA峰值51.4 GFLOPS，源于每个内存访问操作的浮点运算次数较少。这是有限元方法固有的单元级计算的一个众所周知的特点。考虑关键内层循环（第33至41行），它有12次双精度运算（6次加法和6次乘法）以及13次双精度数组项的引用。因此，每字节的计算强度为12/(13*8) = 0.11，与线形图分析报告的FLOP/L1字节值相同。此值随块大小不变。作为参考，该图还显示了在WARP3D执行（单线程）的方程求解过程中，线形图分析显示的MKL BLAS中性能最佳的例程。

基准测试的几个最终要点

• 代码分析的线形图标签列出了：对于第33至41行的循环，FMA是唯一的“特性”，向量化效率为81%（标量计算的3.24倍增益），向量长度= 4个双精度浮点数（AVX2）。此外，内存访问模式的线形图标签确认了例程中两个内层循环的所有数组都采用了步长为1的访问。
• 内存对齐。编译过程使用了–align array64byte选项，该选项对于AVX-512指令的处理器获得最佳性能是必需的。对于此处使用的支持AVX2的处理器，32字节对齐即可。尽管如此，由于驱动程序.f和例程.f文件的单独编译代表了构建过程，编译器无法感知子例程bdbtSoA的正式数组参数的实际对齐。Fortran上下文中的这些可调数组向编译器保证数组元素在内存中是连续的，但不一定具有特定的对齐。注释源代码列表显示了预期的所有数组的未对齐访问。

  对!DIR$ ASSUME_ALIGNED b(1,1,1):64, d(1,1,1):64类型或单独在两个内层循环之前使用!DIR$ VECTOR ALIGNED的源代码指令进行的实验，并未带来可测量的性能提升。编译器列表显示所有数组都已对齐访问。

  解释可能在于Intel Advisor的消息，即所谓的剥离循环和剩余循环（用于执行初始循环迭代直到访问对齐）即使在没有源代码指令的情况下也没有执行。在运行时，由于编译驱动程序.f时使用了–align array64byte选项，进入循环前数组的对齐很容易确定为所有情况下的64字节。

单位步长。一个简单的实验表明了步长为1的数组引用对性能（至少在此处理器上）的关键重要性。在第34行，将数组ek_symm(span,utsz)的设计/使用切换为ek_symm(utsz,span)，并将ek_symm(i,j)更改为ek_symm(j,i)，导致性能从8.08 GFLOPS下降到2.39 GFLOPS，下降了70%。

摘要

向量化为包含大量基于数组的计算的代码提供了潜在的加速——这些加速放大了通过更高层次的并行计算（使用线程和集群上的分布式执行）所获得的性能提升。向量化的关键特性包括可调的数组大小，以反映各种处理器缓存和指令能力，以及内层循环中的单位步长访问。

随着硬件设计者在新兴处理器上将向量寄存器的数量增加到八个双精度浮点数（并希望能更多），向量化对提高性能的重要性将继续增长，以克服性能增长停滞的时钟速率和线程可扩展性问题。

Intel Advisor，特别是新的线形图分析功能，为不同专业水平的代码开发人员提供了一个方便的GUI界面，其中实现了相关的性能信息和建议。

致谢

作者非常感谢他过去在全球范围内任职的研究生、博士后和其他合作者为使WARP3D成为研究人员有用的代码所做的许多重要贡献。有关完整列表和赞助商，请参阅warp3d.net。这些提高代码性能的最新努力是与Kevin Manalo博士和Karen Tomko博士合作完成的，并得到了俄亥俄超级计算中心的Intel®并行计算中心的支持。

参考文献

"Intel Advisor线形图分析," The Parallel Universe, 第27期, 2017。

"使用Intel® AVX-512提高性能的向量化机会," The Parallel Universe, 第27期, 2017。

"Fortran中的显式向量编程，" Corden, Martyn (Intel), 2016年1月。

作者简介

Robert H. Dodds Jr.博士，NAE，是伊利诺伊大学厄巴纳-香槟分校的M.T. Geoffrey土木工程主席（荣誉）。他目前是田纳西大学诺克斯维尔分校土木工程系的特聘教授。

尝试Intel® Parallel Studio XE中的向量化工具

下载30天试用版