朱文强，傅 游，梁建国，郭 强，花 嵘

1(山东科技大学 计算机科学与工程学院，山东 青岛 266590)2(山东省计算中心(国家超级计算济南中心),济南 250101)

E-mail:huarong@sdust.edu.cn

1 引 言

分子动力学(molecular dynamics,MD)[1]作为一个科学计算领域的重要应用，它是高性能计算领域中具有很强代表性和挑战性的问题之一.针对不同的体系结构，研究者在HPC平台上提出各种技术进行分子动力学的加速运算.

2007年，James N.Glosli等[2]利用分子动力学实现了Kelvin-Helmholtz(KH)不稳定性的首次微米级模拟，进展主要体现在容错、内存优化和高效并行3个方向；2009年，David E等[3]针对分子动力学模拟生物分子系统，制造了一台特殊的超级计算机Anton，特殊的硬件极大地提高了MD的计算速度；张勤勇等[4]分析讨论了分子动力学模拟的算法特征和计算特点，采用基于消息传递的MPI设计平台,在可扩展机群上实现了并行化,获得了90%以上的并行效率；赵文辉等[5]提出一种针对分子动力学的Atomreorder算法，通过重新调整粒子顺序来增强程序的局部性，进而减少Cache的丢失.

随着图形处理单元(Graphics Processing Units，GPUs)功能的增强，加之分子动力学大多为数据并行计算，非常适合在其上运行.Jens Glaser等[6]在GPU上对HOOMD-blue程序进行并行优化，最终较CPU取得了12.5倍的加速；杨决宽等[7]在GPU平台提出了一种加速分子动力学仿真计算的方法，并计算了固体氩的导热系数，获得10到11倍的加速；张帅等[8]在GPU平台将分子动力学模拟中的原子划分和空间划分相结合，优化了用于近程力计算的Cell Verlet算法.目前常见的分子动力学软件NAMD,VMD,LAMMPS,AMBER,GROMACS等均已发布支持GPU计算的版本.硅材料作为半导体和太阳能电池行业的重要材料之一，非平衡分子动力学模拟是研究硅晶体热物理性质的有力工具，克服了体状硅和硅纳米线导热系数的有限尺寸效应[9].针对硅晶体分子动力学模拟，中科院过程所的侯超峰等[10,11]在GPUs上实现了优化.

在“神威·太湖之光”[12]上开展分子动力学研究的单位也较多.马跃等[13]提出一种基于收敛策略的国产CPU性能测试模型；中国科学院过程工程研究所在此系统上对Cahn-Hilliard方程进行并行优化，将规模扩展到了1，060万核，双精度浮点运算获得了49.39PFlops性能[14]；中国科学技术大学将纳米级粒度分子动力学软件NAMD在该系统上进行移植和优化，针对从核的结构和主从核执行关系对应用进行优化，优化后单个核组的性能较Intel Xeon E5-2650 v2提高3倍，但在扩展性方面，至多可达到325万核，仍需解决多节点性能与可扩展性问题以及整机优化、调度与负载平衡策略等问题[15]；湖南大学从分子动力学算法出发，针对SW26010的体系结构特征，提出了一种基于簇粒子的新算法，然后逐步实现了3种仿真优化方法：SW26010的并行扩展、内存访问优化和矢量化.经过优化处理后，在单个计算节点上实现了14倍加速，实验中应用了近24万个计算节点(60万核)，获得了近似线性的加速[16]；山东大学在该系统上进行了分子动力学软件LAMMPS模拟，从4个方面进行性能提升，优化后单个计算节点的性能相当于一百多个Intel Xeon E5-2650 v2的性能，使用16,384计算节点获得了2.43Pflops性能[17]；中国科学技术大学在该系统上重构分子动力学应用程序GROMACS，主要也是解决内存带宽限制的难题，最终单核组较主核版性能提升16倍，并且设法实现了300万原子核移到798720个计算节点[18].

以上分子动力学程序的移植与优化都取得了较好的效果，但是本文是针对硅晶体分子动力学模拟程序的移植与优化，原子结构以及模拟算法都有较大不同，另外，硅晶体分子动力学模拟使用了更为复杂的多体Tersoff势[19-21]，使得程序在SW26010处理器上移植时内存受限问题变得更加突出，移植难度加大.徐阳等[22]利用OpenACC将Silicon-Crystal应用移植到神威太湖之光上，并进行了数据流驱动任务图并行化，但他们主要偏重于并行编程模型的研究，而本文主要研究的是Silicon-Crystal应用基于MPI+Athread的并行移植与优化，主要是通过将多种优化方法相结合来克服SW26010内存受限问题.

2 背景介绍

2.1 SW26010处理器

“神威·太湖之光”计算机系统是由国家并行计算机工程技术研究中心自主技术研制的超级计算机，该计算机系统峰值运算速度为125.43TFlop/s，功耗为15.37KW，整机内存总容量为1024TB，访存总带宽为4473.16TB/s.计算结点采用的是片上计算阵列集群和分布式共享储存结合的SW26010异构众核处理器，处理器的频率为1.45GHz，包括4个管理单元MPE、4个计算单元CPE及4个内存控制器(Management Controller，MC)单元组成，其中每个CPE单元由8×8阵列的64核心组成，所以SW26010处理器共260个核心(4*64+4=260)；每个核组有8GB本地内存，4个核组共32GB内存，每个运算核心有64KB的局部存储器LDM，运算核心可以离散访问主存，也可以通过 DMA方式将主存数据批量获取到LDM中，同时运算核心阵列同行或同列还可以采用寄存器通信方式进行通信.

SW26010处理器加速线程库(Athread库)是针对主从加速编程模型所设计的程序加速库.用户通过调用Athread接口函数管理从核的计算和访存，每个线程的创建、回收、调度控制和中断异常管理等操作都需要用户调用接口实现，需要用户更精细地编写代码，以避免造成线程阻塞.该方式提高了执行效率，但也增加了调试难度和出错率.

2.2 Silicon-Crystal应用介绍

Silicon-Crystal应用来源于晶体硅的结晶生长.随着半导体行业以及太阳能电池行业的发展，对晶体硅的数量和质量要求不断提高，研究硅晶体生长以更加精确地指导硅晶体的生产具有重要的意义.硅晶体的每个原子周围有4个近邻的原子，4个近邻原子之间构成正四面体结构，中心原子处在正四面体结构中心，其晶体结构如图1所示，硅晶体由这样的四面体累积起来构成金刚石结构.这4个邻居原子的位置各自处于正四面体的4个顶角上，顶角上的原子与中心原子分别提供一个价电子为这两个原子共同拥有，从而形成共价键，每个硅原子和周围的4个原子组成4个共价键，所以硅晶体属于原子晶体.

图1 硅晶体结构示意图Fig.1 Schematic diagram of silicon crystal structure

Silicon-Crystal应用运用MD方法分析了目前比较常用的势函数在模拟硅晶体生长方面的差异，最终选定了Tersoff势来描述硅晶体原子间的相互作用.Silicon-Crystal应用的主要计算核心为计算Tersoff势的TersoffPotent函数,Tersoff势是由Tersoff根据量子力学中键序的概念，基于简单的量子力学理论而提出的原子间相互作用势模型，为一种超级势，能较好地描述硅的非正面体结构，Tersoff势的主要表达式为：

(1)

(2)

VR(r)=Aexp(-λ1r)

(3)

VA(r)=-Bexp(-λ2r)

(4)

其中，fc表示截断函数，Rij表示截断半径，rij表示原子i和j的距离，VR表示邻居原子间的排斥力，VA表示邻居原子间的吸引力，bij表示原子i与原子j间的键级.其它参数如表1所示.

表1 Si的Tersoff势函数参数表Table 1 Major parameters for Tersoff potential function

3 Silicon-Crystal应用在SW26010上的移植

3.1 Roofline模型

Roofline模型是一个可以将浮点性能、计算密度和存储性能等关联起来可视化的性能分析模型.针对特定硬件平台的Roofline模型的建立可以指导该处理器上不同类型应用程序的优化.

Roofline曲线取决于浮点性能峰值和内存理论带宽，本文针对SW26010处理器的浮点性能峰值和访存带宽分别进行测试和计算.

SW26010主从核均支持256b向量化指令，以及乘加融合指令，每个从核支持一条浮点数双精度流水线，主核支持双浮点数流水线.由此可得，SW26010的双精度理论峰值性能为:

主核：4×1.45×8×2=23.2GFlops

从核：64×1.45×4×2=742.4GFlops

SW26010处理器：4×23.2+4×742.4=3.06GFlops

从以上计算结果可知，SW26010的算力主要集中在从核，所以SW26010上进行应用的移植主要是要发挥其从核的性能，因此接下来主要刻画SW26010从核的访存带宽.从CPE角度看，CPE阵列的内存层次结构包括3个级别：与MPE共享的8GB主存、64KB SPM和寄存器.CPE阵列访存方式主要有两种:1)通过具有全局内存地址的常用装入/存储指令(gld/gst)离散访问主存；2)通过直接内存访问DMA方式访问主存.

通过STREAM Triad基准测试测得访问SPM的峰值带宽最高为154.0GB/s，DMA方式访存峰值带宽为90.4GB/s，MPE访存峰值带宽为39.7GB/s，gld/gst方式访存效果最差，峰值带宽仅为5.92GB/s.

根据上述数据，可以得到SW26010处理器的Roofline性能曲线如图2所示.

3.2 移植原则

结合SW26010处理器的Roofline模型与Silicon-Crystal程序的特征，确定了移植过程中4条需要遵循的原则.

图2 SW26010的Roofline模型Fig.2 Roofline model of SW26010

1)Silicon-Crystal应用是一个计算密集型的应用，从Roofline模型构建过程可知，SW26010处理器的主要算力在从核阵列，因此，应尽量将应用的热点函数放在从核阵列上进行计算，同时保证从核阵列负载平衡.

2)从Roofline模型可以看出,SW26010有限的内存带宽(90.4GB/s)导致较高的每字节浮点运算比率(33.85Flops/Byte)，超过 KNL(7.05Flops/Byte)的4倍多，而其处理器峰值性能(3.06TFLops)却与Intel KNL 7290(3.46TFlops)相当，这意味着SW26010处理器强大的运算能力会受到内存带宽的限制.为避免使应用变成访存受限，必须充分利用内存带宽.为了最大程度地减小内存访问开销，应尽量限制随机数据访问.

3)Silicon-Crystal应用在计算过程中存在数据之间的依赖，为了使程序能够更好的并行执行，应在并行前将依赖消除，同时，为提高数据的利用率和传输效率，尽量将所需数据提前准备好.

4)SW26010处理器从核没有超越函数的处理单元，从核在进行超越函数计算时存在离散访问主存的情况，需要用合理的方式避免从核离散访存导致的计算性能下降.

3.3 主从核移植

将Silicon-Crystal应用移植到SW26010处理器，首先移植到SW26010处理器的单个主核上，然后设计合理的方案将热点函数移植到从核进行并行加速.Silicon-Crystal应用移植工作主要分为3个步骤：主核移植、数据处理和从核移植.

3.3.1 主核移植

Silicon-Crystal应用程序是基于C++开发的，SW26010处理器的主核可以直接对C++代码进行编译，但从核却不支持C++代码的编译，而且主核上的C++编译器sw5CC不支持SIMD库.因此，本文用C语言对Silicon -Crystal程序进行重构.完成主核移植后，使用神威平台自带的gprof工具进行测试，发现TersoffPotent函数和pow函数的总运行时间占程序总时间的91%，是程序的主要热点，接下来将围绕这两个函数进行从核并行化优化.

3.3.2 数据处理

通过对初始数据中原子的位置及编号分析发现原子的编号呈现乱序，考虑到从核移植需要利用DMA方式将主存数据取到从核的LDM中，而利用DMA方式对主存进行访问时，访问的是一块连续的主存地址，所以，需要将乱序的数据进行重新排序，以保证数据的连续性，提高后续从核移植时DMA效率.

Silicon-Crystal应用程序的计算数据是由N个乱序硅晶体原子信息构成，其信息包含原子位置、速度、加速度、临时加速度以及4个邻居的编号，通过原子的位置信息可以得到所有原子均匀的分布在一个长方体空间内，其局部结构如图3(a)所示，每个原子的邻居分布在原子周围，原子与4个邻居之间都构成一个正四面体，结合这种特殊的结构，将数据按照原子位置所在的空间坐标系先x，再y，后z的顺序进行编号，如图3(b)所示，整理得到有序的原子数据，同时从图中可以看出，z方向上位移相同的原子会均匀的分布在一个平面上，且每个平面的原子编号是连续的，因此计算数据就可以看成是由多个这样的平面组成的空间结构.

图3 原子与原子邻居之间的空间结构图Fig.3 Spatial structure between particles and their neighbors

计算数据重排序后的原子位置、速度、加速度、临时加速度信息没有发生改变，但是4个邻居的编号发生了改变，为了防止原子与邻居之间相互关系遭到破坏，用排序后的编号对应更新邻居表中排序前的编号，使每个原子都能通过邻居表找到邻居排序后的编号，所有的排序工作结束后，通过对每个原子的邻居编号分析并结合上一段中由多个平面组成的空间结构发现，每个原子的4个邻居位于原子所在面的前一面和后一面各两个，如图3(b)所示，31号原子的邻居分别为前一面的13、16号和后一面的49、52号，而且相邻编号的原子对应的邻居编号也是相邻的，这是由硅晶体内部原子间的结构决定的，通过排序可以清楚地体现出数据的局部性和连续性.

3.3.3 从核移植

本节将对两个热点函数移植到从核进行并行加速，主要步骤分为任务划分与映射和依赖分析与消除.

1)任务划分与映射

Silicon-Crystal应用在计算时，每个原子的计算量大致相同，进行从核移植时，为了从核之间负载均衡，将原子平均分配到64个从核进行计算，又因为硅晶体原子在空间上有较强的局部性，所以选择采用分块的方式，将原子按空间划分为64块，每个从核计算一块，充分利用原子的空间局部性.

以131072个原子为例，每个从核负责计算131072/64=2048个原子，经过3.3.2节对数据进行处理后，原子根据Z方向上位移不同，划分为128个平面，每一个平面包含1024个原子，任务划分后，每个从核负责计算2048个原子，表明每两个连续平面上的原子为一块，分配到对应的从核进行计算.

2)依赖分析与消除

硅晶体原子之间的作用力是相互的，Tersoff势求的是邻居原子对原子本身的作用力所产生的加速度以及原子本身对邻居原子的反向作用力产生的加速度，最终加速度是将原子所受的加速度以及反向加速度累加，由于不同原子的邻居原子之间存在交集，移植到从核后邻居原子的计算与原子本身的计算可能不在同一个从核上，如果两个从核同时操作一块主存空间，会导致数据出错.例如A原子的邻居是B原子，同时A原子也是B原子的邻居，当原子A与B在分配到两个不同的从核上计算时，两个从核都要对主存上的A、B原子信息进行更改，此时就会产生A、B原子的写后写依赖.为了消除这种数据依赖关系，首先，通过变量重命名法将作用力与反作用力所产生的加速度实现分离，在从核中保存反作用力产生的加速度的副本而不直接更新加速度；其次，在主存开了64份冗余空间存放从核计算的反作用力产生的加速度，在从核计算结束后将数据传回主存，由主核将其累加求得最终加速度，实现主核异步串行数据更新，从而避免写后写数据依赖.

4 性能优化

4.1 LDM优化

Silicon-Crystal应用在进行Tersoff势函数的计算时需要用到许多的计算参数，其中就包含表一中的A(eV)、B(eV)、β、n、c等，如果在从核利用gld/gst直接访问主存中的这些参数，会导致访问的延时较大，而且多个从核同时访问时易出现拥堵现象，不利于程序的加速.而每个从核都有各自的LDM局存，从核访问LDM的延时要远小于从核离散访问主存的延时，所以本文在开启从核前，将从核计算所需要的计算参数以及常量整理打包，集中存放到一个结构体中，从核开启后通过DMA方式将数据一次取到LDM中，从核计算时直接从LDM中获取计算所需参数，减少了大量的离散访存，加快了程序的运行速度.

4.2 DMA优化

在SW26010处理器中，主核MPE与从核簇CPEs共享8GB的存储器，从核可以通过两种方式访问主存：1)利用gld/gst直接零散的访问主存数据；2)利用DMA方式将主存数据批量传输到LDM；DMA方式通常延时较大，但是带宽的利用率较高.针对Silicon-Crystal传输数据量大的特点，应尽量使用DMA方式来提升访存带宽.每个原子的信息封装在一个结构体中，结构体包含位置、速度、加速度、临时加速度、邻居编号等信息组成共16元组(x,y,z,vx,vy,vz,ax,ay,az,tmpax,tmpay,tmpaz,nbrmark1，nbrmark2，nbrmark3,nbrmark4)，计算时需要将计算原子传输到LDM中，一个原子信息就要占用112个字节，64KB的LDM空间最多存放580个原子数据，所以需对计算数据进行分块，取512个原子为一块，通过DMA方式将原子信息传输到LDM，既可以减少离散访存延时，又能将一块连续的数据取到LDM中.

考虑到LDM一次存储的数据量有限，在Tersoff势的核心计算过程中取到LDM中的原子信息只用到了原子的位置以及邻居编号，为了减少多余的数据放进LDM中，在主核进行数据处理时将结构体数组(Array of Structure,AoS)转为数组结构体(Structure of Array,SoA)，把所有原子的位置、速度、加速度、临时减速度以及邻居编号单独存放到对应的数组中，如图4所示，从核可以只获取原子的一个参数数据，减少了大量的不必要数据放进LDM中，由原来的16元组减少到了现在的7元组，从核通过DMA的方式可一次传输1024个原子数据，提高了LDM利用率和DMA效率.

图4 结构体数组(AoS)转数组结构体(SoA)Fig.4 Transform AoS to SoA

4.3 软件cache优化

在实现上述优化策略后，中心原子可以通过DMA方式将原子信息取到LDM中，但是原子的邻居编号需要通过查询邻居表才能获得，无法直接通过DMA方式取到LDM中.但是由3.3.2节可知每个原子的4个邻居分布在原子所在面的前一面和后一面各两个，而且相邻编号的原子对应的邻居编号也是相邻的，即相邻原子的邻居之间存在较强的空间局部性，为了减少离散访问原子邻居带来的长延时，本文通过在从核上设计软件cache来充分利用数据的空间局部性，同时用DMA方式代替离散访存.

软件cache的设计主要考虑数据的预取、替换策略、重用策略3方面，但在SW26010上设计软件cache还需要考虑占用LDM的空间以及符合Silicon-Crystal应用特点，既要保证计算所需数据在LDM中，还要保证DMA的传输效率以及LDM的利用率，下面主要从cache大小的确定、数据预取、数据替换策略以及数据重用策略4方面进行详细设计.

cache大小的确定：经过DMA优化后，从核LDM需要存储1，024个中心原子位置、1，024个中心原子的邻居编号等数据信息，剩余空间不足20KB，所以软件cache大小应该小于20KB.硅晶体结构中每个原子有4个邻居需要分别存放在不同的cacheline中，所以将cacheline数量设置为4，cacheline大小根据原子的空间局部性确定，考虑到数据处理后相邻编号的原子对应的邻居编号也是相邻的，而且呈现每32个邻居编号组成一个连续的数据块，而32个原子占用768字节，所以取的原子数应为32的倍数，为了充分利用LDM空间,cacheline大小确定为3KB，最终cache大小为12KB软件cache的数据预取：Tersoff势的计算过程中将数据的前后、上下、左右平面视为对接，在3个维度上构成回环，使得一个维度上前端原子的邻居反而在末端，所以当数据未命中对数据进行预取时，将首地址设置为未命中原子编号向前取临近32倍数的原子编号地址，通过DMA方式一次取连续的3KB数据放到一个cacheline中.

软件cache的替换策略：因为原子有较强的局部性，取到cacheline中的原子会按序访问，最先取到cacheline中的原子最先访问结束，所以选择采用 FIFO替换策略来保证替换进来的数据得到充分的利用，同时也降低了软件cache设计的复杂度.替换时将cacheline中第一个原子编号存放到对应cacheline的tag中，计算时可以将邻居编号与 cacheline的tag进行比较，快速确定是数据否命中.

软件cache的重用策略：数据整理后已知4个邻居分布在上下两层，而且根据每个维度的大小原子邻居呈现周期性变化，前一个周期所有原子的2号和4号邻居原子会被作为后一个周期所有原子的1号和3号邻居原子再次参与计算，即邻居面上同一个原子会被作为中心面上两个原子的邻居参与计算，为了提高cache中数据的重用率，采用了全相联映射机制，降低cache设计的复杂度.

4.4 从核超越函数优化

完成以上优化后，利用penv_slave2_gld_count接口测得从核存在大量离散访存，主要原因是神威平台的编译环境中默认链接的数学库需要离散访存操作，而程序中存在大量的pow、exp、sin以及cos等数学库函数的使用，导致大量离散访存，减低了性能.如果改用在CPE上实现的已经优化的GNU数学库，虽不会产生gld/gst指令，但要通过查LDM中的相应函数表实现，sin和cos函数表共享4KB空间，而pow和exp的函数表大约占用9KB空间，4个库函数共需要占用大约13KB的LDM空间，使得软件cache和各种数据缓冲能够使用的空间减小.为此，本研究使用多项式近似计算超越函数的值.

ex函数求解方法如下：利用式ex=2xlog2e进行转换；令xlog2e=p+q,则ex=2p·2q,其中p为浮点数的整数部分，q为小数部分且0≤q<1；通过指数函数求解2p；利用式2q=eqln2进行转换；最后利用泰勒公式展开的多项式P(x)进行求解,在P(x)的阶数为11时，|P(x)-ex|<10-14可以满足计算精度的需求.

pow函数的求解方法为：x<0时，y为整数，可以直接求解；x>0时，通过式xy=eyln(x)求解，其中，利用公式：

(5)

求得ln(x)，再利用exp和ln函数求解pow函数.

对于sin函数，则利用sin函数的周期性得到sin(x)=sin(x%2π)，令y=x%2π，为了快速收敛，利用sin函数的对称性，即：

(6)

将y化简为区间[0，π/2)内的值，然后利用泰勒展开式：

(7)

进行计算，最终求得sin(x)的值.cos函数计算与此类似.

从核直接调用以上实现的超越函数进行求解，不仅消除了从核调用超越函数产生的离散访存，而且节省LDM空间.

4.5 任务分步流水优化

经过使用软件cache后，极大地提高了数据的重用率，减少了DMA的访存次数，但是由于LDM仅有64KB，依然无法将整面原子计算所需要的数据以及计算后产生的数据完全存放在cache中，而且仍然会有许多原子作为不同原子的邻居多次传输到LDM中，导致传输效率低，所以本文设计了一种分步流水的方式将Tersoff势的计算分为两部分进行，对计算所需以及计算产生的数据分到不同LDM中进行存储.

深入分析Tersoff势的计算过程可以看出在每个原子的计算主要分为两部分，第1部分为截断函数的求解，通过原子与邻居的位置求得相互之间的距离，然后代入2.2节中公式(2)进行求解；第2部分为硅晶体势能求解，将求得的原子与邻居之间的距离代入公式(3)和公式(4)求出原子之间的吸引力和排斥力，结合第1部分截断函数的求解，通过公式(1)得到硅晶体的势能.计算过程中，主要的存储数据为第1部分求解所需原子的位置以及第2部分利用原子之间相互作用力求得的原子的加速度，如果将一整面计算所需要的位置数据存放到LDM中，计算产生的原子加速度将无法存放，因此为了使计算数据能够完整地存储到LDM中，将前后两部分计算分配到两个从核进行计算，中间交互使用寄存器通信，使用从核寄存器的点对点通信全片聚合带宽高达2043.4GB/S，合理使用寄存通信可以进一步减少从核冗余的访存操作，极大提升计算性能.

图5 从核间双流水模式Fig.5 Step-to-step pipeline and double buffering between CPEs

将从核阵列分为偶数核和奇数核，偶数核从主存获取计算所需的原子位置信息存储到LDM中，然后进行截断函数的计算，计算结束后通过寄存器通信传输到奇数核，奇数核拿到计算数据后进行吸引力、排斥力和加速度的计算，将计算后数据存储到奇数核LDM中，在奇数核计算的同时，偶数核可以进行下一个原子截断函数的计算，这样便形成了从核间任务的分步流水计算，如图5(a)所示，此时将计算原子以及邻居原子一次存到偶数核LDM中，将计算后的加速度数据存放到奇数核LDM中，减少DMA次数，而且提高DMA效率.

实现任务分步流水优化后，从核LDM仍然有剩余空间，考虑到每个从核要对多个连续的原子面进行计算，所以使用双缓冲优化将计算与通信进行隐藏，如图5(b)所示，偶数核在计算前发起下次所需数据的DMA请求，到下次计算之前检测数据是否传输完成，奇数核在发起写回主存的DMA请求的同时直接进行下一次的计算，直到下一次计算结束后判断传输有没有结束，奇数核和偶数核都实现了双缓冲流水，结合任务分步流水，形成了从核间分步流水从核内部双缓冲流水的双流水模式

5 性能实验结果与分析

为了测试验证上述研究内容的效果，将Silicon-Crystal应用分别在SW26010异构众核处理器和Intel Xeon E5-2620 v4处理器上的性能进行测试对比，表2列举了试验平台信息.

表2 实验用软、硬件环境配置信息Table 2 Experimental software and hardware environment configuration information

实验将运行在1个主核上的Silicon-Crystal应用作为测试基准，对第4节中的5步优化逐步进行测试.测试数据规模为：131072个原子构成的回环，迭代计算次数为1000次.

5.1 单核组性能优化评估

本文所有优化方法的性能测试结果如图6所示，从图6中可以看出：

图6 多种方法优化后程序的运行时间Fig.6 Running time of the program after implementing each optimization

1)进行简单从核并行后，性能不但没有提升反而有所下降，主要原因是从核的离散访存占时较大，导致计算优势难以发挥.

2)对从核访存方面进行的优化效果进行评估.首先通过数据结构的转换使其更适合DMA传输，提高从核的访存带宽；针对应用中存在的不规则数据访问，结合硅晶体的结构特征进行数据预处理，使数据更具空间局部性，从而利用软件cache进一步提升从核的访存带宽；针对应用中存在的大量超越函数的离散访存，通过在从核上定制的数学库函数解决了其离散访存问题；为了进一步减少从核访存，还通过邻居预取以及任务分步流水方法进一步减少从核访存.最终单个核组的运行时间由778.55s减少到60.37s，其性能比测试基提升了12.89倍.同时还可以看出,SW26010处理器的单核组优化性能相较于Intel Xeon E5-2620 v4处理器串行程序提升了8.7倍.

5.2 多节点性能评估

单核组性能优化及评估结束后，本文又在神威平台上利用MPI+Athread实现了Silicon-Crystal应用的多节点并行版本，并对其进行了强、弱可扩展性测试及分析.

强可扩展性：对Silicon-Crystal应用的数据进行扩展，将数据量扩展到67108864数据，进程的规模最大扩展到了512个核组，测试结果如图7所示.

图7 Silicon-Crystal应用在从核上的强可扩展性分析Fig.7 Strong scalability of Silicon-Crystal application

从图7可以发现随着核数的增长，运行时间在逐渐减少，到最后趋向平稳，因为随着核组数的增长，数据的切分越来越多，进程之间的通信也随之增多，进程间的通信占比越来越大，程序运行时间逐渐受到通信的限制趋向于平稳.

表4 弱可扩展性测试规模Table 4 Scale of weak scalability test

弱可扩展性：随着核组数的增长，数据也跟着同比增长，由最开始的4个核组运行524288数据扩展到512个核组运行67108864数据，数据规模以及不同规模下程序的执行时间如表4所示，可以看出随着核组数和数据量的同比增长，执行时间并没有随着通信量的增加而加长，呈现相对平稳状态，表明Silicon-Crystal应用程序具有较好的弱可扩展性.

6 总 结

本文基于国产异构众核处理器SW26010的体系结构特点，使用神威平台的MPI+Athread并行编程模型对Silicon-Crystal应用进行了并行化.在对应用进行主核移植后，针对Silicon-Crystal应用程序热点函数特点，从5个方面进行单核组优化，使其性能相比主核串行基准测试程序提升12.89倍.与Intel CPU平台串行程序相比，提升8.7倍.

下一步将针对多节点运行进行更深入的优化，并对大规模可扩展性进行深入研究.