计算流体力学程序单核指令级优化方法

刘闯,何峰,肖兮,董小社,张兴军

(西安交通大学计算机科学与技术系,710049,西安)

通过计算流体力学程序模拟真实情景是提高流体机械的制造精度和节约成本的重要手段,计算流体力学程序具有计算密集、处理数据量大的特性,因此常需要大量的计算资源。尽管计算机的计算能力有了飞速的发展,计算流体力学程序对计算性能的利用率仍然存在不充分的现象,导致计算资源的浪费。

目前对于流体计算程序性能提升的研究主要体现在并行化或异构加速上[1-3],例如采用多核系统并行处理或采用协处理器加速计算。超级计算机上部署的大型计算流体力学程序大多都采用了上面两种计算方式,但是作为一个重要的影响因素,程序的单核执行效率并没有得到足够的重视,很多程序的单核性能只能达到系统理论峰值的5%～10%[4],导致大量的计算资源被浪费,因此对程序单核性能的优化应引起重视。

针对流体计算程序的优化,一些文献提出了综述性的优化方案和思路[5-7],或针对具体程序进行了优化[8],但是这种优化主要集中于对云计算、网格计算等并行化计算的优化,对单核级优化的工作不够深入;针对单核级的程序优化的研究中,主要优化方案是流水优化或访存优化,优化步骤一般是首先分析程序特性或测试程序热点以找到瓶颈,然后通过调整代码结构流水性能[9]或调整数据的组织方式优化访存性能[4,10-12];从处理器的角度,贺爱香、Oyarzun等根据ARM处理器的计算特性对问题提出优化方案[13-14],这些研究具有一定的通用价值,但是对于大型计算程序常用的Intel处理器并没有深入研究;申小伟等通过对编译器等运行环境的优化来提升程序性能[15],但是这种方法对于跨专业的程序优化人员具有一定的门槛,且在公用的计算环境中难以部署。根据上述背景,本文针对轴流压气机模拟的流体计算程序,提出一种针对计算流体力学程序单核的优化方法。总体优化思路有3个方面:①优化指令结构,减少由相关引起的指令流水断流,使程序充分发挥CPU的计算能力;②调整存储结构和访存结构并进行容器优化,从而增强存储和访存的性能;③减少冗余指令,优化底层代码数量和结构从而减少程序执行所需的时间空间。本文在轴流压气机模拟程序上实现优化方法,在TIANHE-1A超级计算机和商用服务器上进行测试实验,实验结果表明,本文提出的优化方法可以有效地提高程序的缓存命中率、优化流水连续性、减少目标代码数,从而大幅提升程序单核执行效率。

1 流体计算程序简介

1.1 流体计算程序背景

本文采用轴流压气机模拟程序[16-17]对轴流压气机转子进行模拟,程序采用有限体积法进行空间离散,3步Runge-Kutta法进行时间推进,控制方程组采用笛卡尔坐标系下的相对流动控制方程组,形式如下

(1)

湍流模型采用SA湍流模型,为加速收敛程序采用三层网格,通过在不同网格层上轮流迭代从而使不同频率的误差分量得到有效地衰减。

程序模拟36叶片压缩机工作时的流体状态,采用36个进程,每个进程对一个叶片进行模拟,单叶片网格数约为42万,总计约1 535万。在边界处理时需进行通信,通信采用MPI完成。程序主要开销是在单个进程的计算上,通信等待时间相对整个程序时间占比较少,且优化不涉及通信时间的优化,因此并行开销可忽略不计。

1.2 流体计算程序实现细节及性能分析

考虑程序的缓存命中率对程序的存储和计算特性进行分析。存储方面,程序中网格的存储采用STL库中的Vector容器进行存储,采用该容器的主要原因是Vector容器是顺序存储容器,有较好的随机存取性能,以及可以利用STL中提供的丰富的操作,功能强大。根据物理网格模型,容器的维度为四维,存储维度顺序为x、y、z、n,其中x、y、z为空间维度,n为多层网格的层号。

计算方面,程序计算主体的运行时间约占整个运行时间的95%以上,剩余时间做程序流程控制、数据的通信与保存等操作,因此程序优化的主体在计算部分。在计算部分中,程序对整个网格层采用for循环进行更新,大多数循环深度为3层,由于边界上存在虚拟网格,循环次数约为43.5万次,略多于网格数,每次迭代中更新一个或多个变量。此外,为实现复杂计算,程序还调用了math库中的部分计算函数,考虑到计算函数计算开销较大,在某些调用中可采用等价短指令替换。

为进一步分析程序计算性能,本文采用性能分析软件对程序进行分析。结果显示,程序的流水线性能未充分发挥,双向中断率均超过50%,暂停周期数较高,导致执行某个程序的指令平均时钟周期(CPI)与理论值有差距。因此,在优化时应考虑流水优化;缓存命中率约为50%,访存上也有优化空间;此外程序执行指令数为38 534×109条,在优化中需寻找潜在的指令优化来减少冗余代码。

2 优化策略研究

本文所述方法应用在TIANHE-1A超级计算机系统中Intel Xeon X5670 CPU上,采用Intel CPU作为示例是考虑到大多数的高性能计算环境都采用Intel处理器,且Intel CPU中采用的加速部件在其他处理器中有类似部件设计,因此以Intel CPU为基础设计的优化方案具有较好的普适性。

图1 Intel CPU架构示意图

图1给出了Intel CPU的架构示意图,该示意图基于Intel Xeon X5670 CPU所采用的Nehalem架构,主要体现CPU针对单核程序的加速方式,对其余部件做了省略和简化。图1中CPU包含两个核,每个核内设置独立的计算部件,从逻辑上分为取值部件、解码部件、顺序回退部件、计算单元,图1中的箭头代表了部件间的数据通路。基于上述架构,Nehalem CPU架构部署了取指(IF)、译码(ID)、执行(EXE)、写回(WB)四阶段流水机制。此外,Nehalem架构还采用了超标量技术,在CPU中建立了4条流水线,因此理论上可以同时执行4条流水指令。每个核设置独立的L1和L2缓存,多个核间共享L3缓存,方便进行数据缓存,提升访存性能。其中一级缓存分为数据缓存和指令缓存,每个缓存容量为32 KB,采用8×64组相连的映射方式,访存速度为4个时钟周期;L2缓存容量为256 KB,采用8×512组相连映射方式,访存速度为10个时钟周期;L3缓存容量为8 MB,采用16×8 192组相连映射方式,访存速度为40～75个时钟周期。

基于对Nehalem架构的分析,程序优化的思路主要集中在使程序充分发挥流水性能和访存性能上。

除针对CPU架构的优化外,还需考虑程序本身的简洁性。对于不同形式的程序,通过编译器编译得到的目标代码有很大差距,若有冗余和过多同义代码会导致性能下降。为保证程序性能,在编码时需考虑高级语言被编译后得到的目标代码是否简洁高效。

结合上述分析,本文提出的流体计算程序单核优化方案主要考虑3个方面:①优化存储结构和访存结构,发挥高速缓存效率;②调整程序的指令结构,发挥流水性能;③优化程序指令,使程序底层指令数得到削减和优化。前两个方面面向CPU对单核指令执行的高速访存和流水优化部件展开,目的是充分发挥CPU加速部件的加速效果。第3个方面针对程序代码本身优化,考虑到不同指令在底层执行时采用不同的汇编指令翻译,程序的目标代码应足够精简高效。

2.1 访存优化

访存优化的目的是使程序在数据访问时缓存的利用率更高,主要优化手段是调整数据的存储结构和访存结构,首先考虑存储结构。

(a)访存在内存上连续

(b)访存在内存上离散

程序的访存方式是通过循环遍历的方式对网格进行访问,优先访问网格层号n,其次为x、y、z轴,这种存储结构和访存方式的不匹配造成了缓存效率较低。图2给出了不同存储方式下访存情况的对比,对于M×N矩阵,规定访存顺序为从(0,1)到(0,N),若矩阵采取aM,N的形式存储,则矩阵在内存中的物理分布如图2a所示,在访存时数据分布在连续的内存中,缓存会将连续内存存入,在下次访存时可直接取用缓存,减少内存访问。若矩阵采取aN,M的形式存储,则矩阵在内存中分布如图2b所示,在访存时会产生跳跃,因此缓存无法命中,需要进行内存访问,影响程序运行效率。

为解决访存问题,将数据的存储结构调整为an,z,y,x,这样调整的好处有:①存储结构更改后x网格内相邻各点在内存中存储的更近,更有利于程序的连续访问;②对于多维容器来说,由于x>y>z>n,更改后的存储形式减少了容器在较低维度所创建的数量,减少了目标代码的指令数。例如,对于容器a20,10,需先创建一个容量为20的容器,每个容器中再创建一个容量为10的容器,总共创建21个,若改为a10,20,则仅需创建11个容器。

下面根据存储结构调整访存结构,访存顺序的控制由for循环语句进行。为使访存顺序更适合存储结构,将循环结构和存储维度调整为一一对应的关系。即对于an,z,y,x形式的数据,处理循环为最外层循环处理n、次外层处理z、内层处理y、最内层处理x。对不满足这种处理方式的操作,考察操作的相关性,若无相关性可通过拆分循环来优化访存顺序。

访存方式的优化和存储结构的优化有类似的优点:调整访存结构有利于高速缓存;调整后循环的控制变量初始化次数减少,消减底层指令数。

2.2 流水优化

考虑影响流水线性能的问题,在程序一级影响因素主要有:由分支转移引起的控制相关会导致流水失效;由写后读引起的数据相关,这些相关会导致流水暂停,使流水性能下降。Nehalem CPU架构部署了四阶段流水机制,虽然该架构加入了分支预测和数据旁路的方式减少流水线中断,但是仍不足以完全消除相关性引起的性能损失,因此在优化时仍可消除相关性。

考虑计算中的循环部分,由于循环在每次迭代时需进行一次判断,若循环内指令较少会导致计算和分支的比例过低,产生控制相关导致性能下降,同时循环中对同一变量的反复读写可能引起数据相关。针对控制相关和数据相关问题,一般采用循环展开并引入多中间变量的方法进行优化,循环展开可以在以下方面得到优化:①减少分支数;②减少循环控制的指令开销;③减少循环内计算的数据相关。循环展开示例如下

for(i=0;i

s+=a[i];

改为:

for(i=0;i

s1+=a[i];

s2+=a[i+1];

}

s=s1+s2;

展开因子k=2,即每次循环计算次数为两次。两种代码的流水情况如图3所示,流水性能得到了提升。

(a)循环拆分前

(b)循环拆分后

理论上,在展开足够多循环的情况下可以完全消除流水中断,展开因子应大于硬件标量数与执行周期的乘积,如对于具有3个浮点乘器件,每个浮点乘需5个时钟周期的指令集,展开因子应大于15;但是实际操作中需考虑寄存器溢出,即中间变量过多导致寄存器不够,触发访存。为此还需保证最优循环因子k满足

max{k};F(k)≤r,k∈Z

(3)

式中F(k)为k次展开需要的寄存器数。

此外,解除数据相关还可以采用长计算拆分为短计算、创建临时变量的方式解除相关、多个连续的写后读相关语句交叉等方法。

对if语句引起的控制相关,分支预测部件往往不能正确预测走向,预测错误的惩罚往往达到十几个到几十个时钟周期,带来流水性能损失。在X86指令集中,条件数据传输指令可以被编译为普通指令流水线的一部分,在判断失败时没有惩罚,因此采用条件数据传送代替条件控制转移能提升程序性能。常用的条件数据传输指令是三目运算符,条件传输形式如s=a

x=a

x′=b

boolT=(a

if (!T)x=x′

可以看出,这种编译方式将分支转移操作放在伪代码最后一行的if语句中,在目标代码中这个if语句通过条件拷贝指令cmov实现,分支的结果只在这一条指令中体现,与上下文无关,从而消除了分支开销。通过这种转换可以提升分支在流水线中的稳定性,减少预测错误的惩罚。

2.3 指令优化

指令优化的目的是从指令一级减少冗余代码。冗余代码主要指程序编译后的目标代码中产生的冗余,一般主要产生于如下两种情况。

(1)完成相同功能的不同语句在经过编译器编译后会产生不同的目标代码,若选择不合适的语句会引起冗余的目标代码。这种冗余包括两种,一种是目标代码量的增加,指对于同样功能的代码由于实现方法不同导致编译后代码条数不同;另一种是代码复杂度的增加,CISC指令集中的指令执行时间不同,因此同样的指令条数下指令执行时间有所不同。

(2)为实现方法调用,操作系统需保存当前运行时信息,并在运行时栈开辟新的空间,在调用结束时恢复信息并释放空间。这种调度工作会体现在目标代码中,虽然一般认为正常的调用过程对程序性能影响不大,但当调用在程序中占比过大时会影响程序性能。调用引起代码冗余的一种示例如下。

非调用版指令如下

s=a+b;

目标代码为

1 mov %rdi,%rbp

2 mov %rax,rbp

3 add %rbx,%rax

调用版指令如下

Sum (inta,intb)

{returna+b;}

目标代码为

1 push %rbp

2 push %rbx

3 sub %8,%rsp

4 mov %rdi,%rbp

5 mov %rax,rbp

6 add %rbx,%rax

7 add %8,%rsp

8 pop %rbx

9 pop %rbp

10 ret

调用版的第1、2行用来实现现场保存,第3行为新调用栈分配空间,然后执行加法功能,第7行释放调用栈,第8、9行恢复现场,最后返回。可以看出,当实现同样加法功能的s=a+b被替换为方法调用Sum(a,b)后,运行时栈调度的开销超过了方法本身。虽然编译器中有对高级语言的自动优化功能或调用产生冗余代码的问题,如自动将短调用进行内联处理,但是编译器优化的前提是保证程序的安全,因此面对具有内存别名使用或变量类型未知等问题的代码不能进行自动优化,需手动对高级语言进行优化。

分析本程序得出,导致性能下降的主要问题是:①部分库函数调用时产生的隐藏代码增加了代码量,增加包括调用栈引起的代码开销和库函数本身复杂性引起的代码开销;②代码及代码结构本身的计算复杂度较高,因为有较多耗时的除法计算需求。为解决上述问题,考虑使用简单指令代替复杂指令,并结合访存优化及流水优化考虑汇编指令在系统中的运行状态进行优化。

下面讨论指令优化的具体实现。针对库函数调用的问题,结合程序分析发现,程序主要调用的库函数为数学库和容器类库。数学库的乘方、开方等操作一般不可替代,但是对于低阶乘方,可以直接使用数乘代替。如pow(a,2)可以改写为a*=a,这种替换可以减少调用开销,同时简化方法本身的复杂度。

对容器类库函数的修改可考虑对C++语言自带容器数组进行二次封装代替库容器,在本文的程序中为采用数组代替Vector容器。从可行性分析,Vector和数组有较为类似的存储结构,因此在存取功能上具有可替代性。Vector相对数组的优势在于灵活的操作方式和内存分配方式,程序中对容器的主要操作方式为随机访问,此外有少量的动态内存分配操作,因此在二次封装数组时需要为数组增加随机访存函数,从而实现对Vector的替代,所以用数组代替Vector具有可行性。从必要性分析,Vector是STL封装的容器库,在使用时涉及库调用会增加底层代码量,容器的存储结构上,Vector为实现动态增长,在分配空间时会分配额外的内存,在多维空间上,这种存储空间的浪费变得更加明显,在访存时也会增加额外开销,高速缓存的命中率也会下降。

针对代码结构结合访存特性考虑编译器对语句的处理方式,由于对变量的访问常被编译为寄存器访问操作,对引用的访问被编译为内存访问操作,因此为优化访存性能,需将循环中频繁引用访问修改为变量访问。针对代码本身,考虑科学计算中常出现的乘除法指令,在Intel提供的参考机中部分算数指令性能见表1,其中延迟表明实际运算周期数,发射表明两次运算最小间隔周期数,容量表明该功能的运算单元数。

表1 Intel参考机算数指令性能

由表1可以看出,乘除法运算执行时间较长,除法尤为明显,因此对乘除法的优化也应加以考虑。针对非浮点数的乘除法指令,采用移位运算代替乘除法的方法进行处理。该方法在乘除法数量较多且乘数或除数较为规律时可产生较好的优化效果。对于浮点数除法,可采用公共除法移出的方法将除法转换为乘法,减少执行时间。

3 实验结果

对第一节中的流体计算程序采用上述优化方案进行优化,并在商用服务器对各版本程序采用perf性能检测工具对程序进行小规模性能分析。商用服务器采用Intel Xeon CPU E7-8850 @ 2.00 GHz(8×10核)CPU,内存125 GB,编译器采用gcc 4.8.5。在三层网格上迭代轮次为5、5、10。程序优化前后的主要指标的分析结果和执行时间见表2。由于对容器以及对指令的优化,程序所使用的运行时空间的大小也有所下降,内存占用实验结果见表3。

表2 程序主要指标优化前后的比较

表3 优化前后占用内存情况

由表2、3可以看出,程序性能的主要提升在于优化后导致指令数的减少,在系统级优化后缓存性能和流水性能也有所提升。在规模为5、5、10的程序中,经过优化的程序每周期指令数提升24.24%,指令暂停周期下降25.86%,流水中断率下降6%,缓存命中率提升14%,执行时间下降70.56%,空间消耗下降49.97%。

将程序部署在TIANHE-1A超级计算机系统中测试性能,采用4节点运行36进程。采用的CPU为Intel Xeon X5670(2×6核),内存24 GB,编译器为icc_xe_2013.0.07。三层网格迭代轮次采用50、50、10 000作为程序完整规模运行的实验,测试结果见表4。

由表4结果可以看出,在TIANHE-1A系统上,三层网格迭代轮次分别为50、50、10 000计算规模时,程序耗时下降68.34%,空间消耗减少55.43%,与商用服务器上测试结果类似。根据测试结果可以看出,本文的优化方法能有效提升程序对CPU计算能力的发挥,大幅减少程序的底层代码,在测试程序中取得较好的优化效果。

表4 程序完整规模运行的实验测试结果

4 结 论

根据计算机系统单核级指令执行的加速特性和编译特性,提出一种针对计算流体力学程序的单核优化方案,主要优化方面包括:①优化指令结构,减少指令流水的断流,使程序充分发挥CPU的计算能力;②调整存储结构和访存结构,增强存储和访存的性能;③减少冗余指令,优化指令结构,从而减少程序执行所需的时间开销和空间消耗。

将本文方法应用在压气机转子模拟程序中,测试结果显示优化方案能有效减少程序指令数、增加缓存命中率,减少指令流水断流,在TIANHE-1A计算机系统中进行测试,结果显示优化后时间性能提升约68.34%,空间性能提升约55.43%,加速比提升3倍以上,具有较好的优化效果。