CN112100099B - 一种面向多核向量处理器的格子玻尔兹曼优化方法 - Google Patents

Abstract

本发明公开了一种面向多核向量处理器的格子玻尔兹曼优化方法，包括以下步骤：S100、通过设定的预取、计算和写回三个状态针对芯片内核的向量访存部件进行数据传输时间隐藏的设计；S200、设计基于表面切割的复杂标量程序的向量化方法，以及面向基于向量的快速指数算法，并使用改进的三维格子玻尔兹曼方法进行验证实例验证，实现效率提升；S300、基于单核的验证流程自适应到多核，实现快速并行计算。相比现有技术，很好地解决了向量的并行运算的问题，缓解了DMA带宽低导致的时延长，实现了指数运算的高度向量并行计算。

Description

一种面向多核向量处理器的格子玻尔兹曼优化方法

技术领域

本发明属于高性能计算领域，具体涉及一种面向多核向量处理器的格子玻尔兹曼优化方法。

背景技术

围绕高性能自主加速器关键技术展开研究是各国超级计算的热潮。多核向量处理器高性能计算的难点主要体现在：设计多核向量的并行算法难；数据量大且DDR访存带宽低延时长；exp等库函数缺乏高度并行的向量化版本。在算法程序进行向量计算时，需要先把数据从DDR读入到核内向量地址空间(AM)，才能进行计算。DDR的存储空间大，但读写速度明显慢于AM，降低了运行效率。

格子玻尔兹曼方法(LBM，Lattice Boltzmann Method)是一种应用于计算流体动力学的数值仿真方法。它基于量子统计力学的理论，将流体看做大量有质量无体积的微粒，可以通过统计这些微粒的碰撞规律来描述流体的真实流动。格子玻尔兹曼方法主要由格子(离散速度模型)、平衡态分布函数、分布函数的演化方程组成，采用不同的格子、平衡态分布函数以及演化方程，一般便能解决各式各样的流体问题。LBM适合并行运算，但计算机求解流体流动过程中LBM用到取模操作，在转化成向量运算时，取模操作会使得读取数据不连续，难以并行化处理，降低了计算效率。求指数(exp)的原理是通常利用泰勒展开式，存在大量的除法运算，同时，当变量为负数时，泰勒公式中的n需要取很大的值，才能获得较精确的结果，两者都延长了计算时间。

发明内容

本发明提供一种面向多核向量处理器的格子玻尔兹曼优化方法，用以解决现有技术中向量并行算法效率低，数据量大，DDR访存时间长，exp等向量库函数并行度低，导致计算效率低等问题。

本发明的目的通过以下技术方案实现：提供一种面向多核向量处理器的格子玻尔兹曼优化方法，包括以下步骤：

S100、通过设定的预取、计算和写回三个状态针对芯片内核的向量访存部件进行数据传输时间隐藏的设计；

S200、设计基于表面切割的复杂标量程序的向量化方法，以及面向基于向量的快速指数算法，并使用改进的三维格子玻尔兹曼方法进行验证实例验证，实现效率提升；

S300、基于单核的验证流程自适应到多核，实现快速并行计算。

作为进一步的改进，所述步骤S100具体表现为：根据设定的预取、计算和写回三个状态依次轮换，将DDR数据读写时间隐藏在DSP计算过程，实现对芯片内核的向量访存部件的访存延迟隐藏。

作为进一步的改进，所述步骤S200中使用改进的三维格子玻尔兹曼方法进行验证实例具体按照体-面-边-点的维度分别进行向量计算，从而消除标量的取模操作。

作为进一步的改进，所述按照体-面-边-点的维度分别进行向量计算，从而消除标量的取模操作过程分解为如下步骤：

S201A、计算腔体内部，即计算非固壁网格点；

S202A、计算6个面的内部：先计算6个面内部的网格点，再计算网格点的数据碰撞和迁移；

S202A、计算6个面的交线：通过VLR向量长度寄存器相关函数，设置好需要计算的向量单元长度，再进行各交线上网格点的数据碰撞和迁移；

S204A、计算各条线的交点：设置VLR长度为1，然后通过VLR向量长度寄存器相关函数，设置好需要计算的向量单元长度，再进行各交点上网格点的数据碰撞和迁移。

作为进一步的改进，所述步骤S200中设计面向基于向量的快速指数算法的具体过程如下：

S201B、针对求指数操作，将标准的泰勒展开式，改为n个连续的乘加操作；

S202B、设置多级系数字典，计算时通过查表取值；

S203B、利用指令延迟槽的循环展开。

作为进一步的改进，所述步骤S202A中再计算网格点的数据碰撞和迁移具体表现为：

先计算每个网格的下一个时刻的互作用力时，通过缓存中间结果减少内存和外存之间的数据交换；

再复用网格的当前状态和下一时刻状态的存储空间，实现矩阵复用设计。

作为进一步的改进，所述步骤S300基于单核的验证流程自适应到多核，实现快速并行计算的实现过程为：根据任务的大小和芯片核的个数，动态申请、自动分配DSP的数量，假设芯片核的个数为M，步骤如下：

Step1：根据芯片核的个数M，分别在x、y、z维度上分割为mx、my、mz个子腔体，满足mx×my×mz＝M；

Step2：计算a＝M/mx，b＝M/my，c＝M/mz，a,b,c分别为x、y、z维度上区间的长度；

Step3：针对x维度上长度为M的范围，分成了[0,a-1],[a,2a-1],......,[(mx-1)a,mx(a-1)]的mx个区间；每一个x区间对应的y维度上可划分[0,b-1],[b,2b-1],......,[(my-1)b,my(b-1)]的my个区间，每一个y区间对应z维度上可划分[0,c-1],[c,2c-1],......,[(mz-1)c,mz(c-1)]的mz个区间；

Step4：根据Step3的区间划分将总任务分配到M个芯片核中进行计算。

本发明面向多核向量处理器的格子玻尔兹曼优化方法，首先，通过预取、计算和写回三个状态的设计实现了芯片内核的向量访存部件(即DMA，Direct Memory Access，直接存储器访问)数据传输时间的隐藏；其次，设计基于表面切割的复杂标量程序的向量化方法，以及面向基于向量的快速指数算法，并使用改进的三维格子玻尔兹曼方法进行验证实例验证，修改内层循环，消除标量间的取模操作，实现效率提升，其中矩阵复用设计很好地提升了计算效率，尤其，针对指数运算执行效率不高的问题进行数学形式的转变，并给出多级系数字典；最后，基于单核的验证流程实现自适应多核运算。相比现有技术，本发明针对多核向量程序的难点，解决了向量的并行运算的问题，缓解了DMA带宽低导致的时延长，实现了指数运算的高度向量并行计算。

附图说明

利用附图对本发明作进一步说明，但附图中的实施例不构成对本发明的任何限制，对于本领域的普通技术人员，在不付出创造性劳动的前提下，还可以根据以下附图获得其它的附图。

图1是本发明实施例中芯片的总体框图。

图2是本发明实施例中单DSP内核结构框图。

图3是面向多核向量处理器的格子玻尔兹曼优化方法的流程图。

图4是本发明实施例中AM的内部构造示意图。

图5是本发明验证实例的格子玻尔兹曼B3Q19的模型图。

图6是原始的LBM算法流程图。

图7是本发明实施例中改进的LBM算法流程图。

图8是本发明实施例中延迟槽循环展开示意图。

图9是本发明实施例矩阵复用的设计原理图。

图10是本发明实施例中多核计算实现的方法实例图。

具体实施方式

为了使本领域的技术人员更好地理解本发明的技术方案，下面结合附图和具体实施例对本发明作进一步详细的描述，需要说明的是，在不冲突的情况下，本申请的实施例及实施例中的特征可以相互组合。

本发明的验证方法的硬件平台主要由主芯片、存储(DDR、flash)、外设接口(PCIe、I2C、gpio、gMac)、状态监控(温度传感器、jtag)、时钟、电源、GSM(全局共享内存)和DSP等部件组成，其中主芯片是核心，图1包含了主芯片的总体框架，图2为单个DSP(Digital SignalProcessing，数字信号处理)内核结构框架图。每个DSP内核有1个标量计算单元(SPU)、1个向量计算单元(VPU)和1个DMA控制器，其中向量计算单元有3个向量处理部件(M1～M3)，VPU支持SIMD指令。

如图3所示，本发明实施例提供的一种面向多核向量处理器的格子玻尔兹曼优化方法，包括以下步骤：

S100、通过设定的预取、计算和写回三个状态针对芯片内核的向量访存部件进行数据传输时间隐藏的设计；

具体地，在步骤S100的过程中，程序进行向量计算时，需要通过DMA先把数据从DDR读到核内向量地址空间AM再进行计算，但DDR的读写速度慢于AM。AM空间由16个BANK组成，每个BANK又分为4个SRAM。每个BANK可以独立进行读写，一个BANK内部的4个SRAM可以并行，由于每个核只有1个DMA，且每个DMA只有1个读通道和1个写通道，所以DMA不能同时读AM或者同时写操作AM。本发明将AM向量空间每块SRAM设有3种状态，“预取”状态表示DMA正从DDR预取数据到本块AM空间，“计算”状态表示VPU正在使用本块AM空间，“写回”状态表示DMA正将计算结果写回DDR，写回的地址同读取的地址相同，这样能最大限度地隐藏DDR访存延迟，提高程序效率。如图4所示，BANK内部的4个SRAM必须按照预取、计算、写回的顺序，初始化时，当SRAM0执行预取时，SRAM1～SRAM3部分处于等待状态，每个BANK的同一编号的SRAM执行相同的操作，并都满足上述要求。即本发明根据设定的预取、计算和写回三个状态依次轮换，将DDR数据读写时间隐藏在DSP计算过程，实现对芯片内核的向量访存部件的访存延迟隐藏。

需要说明的是，本发明还可以实现在传输过程执行其他指令代码，当需要用到传输的数据时，再通过检测传输完成，取传输的数据。本发明将DMA传输过程的时间给隐蔽了。同时，这样设计预取和写回可以并行，比乒乓缓冲区方案并发度高，适合I/O吞吐量大的情况。

S200、设计基于表面切割的复杂标量程序的向量化方法，以及面向基于向量的快速指数算法，并使用改进的三维格子玻尔兹曼方法进行验证实例验证，实现效率提升；

在该步骤中，使用改进的三维格子玻尔兹曼方法的验证实例，即针对面向向量的多核LBM算法进行了相应的改进。选择LBM算法作为芯片验证算法的主体，是因其网格在3D空间中，各自独立，适合并行和向量化实现。LBM的B3Q19模型如图5所示，某时刻的平衡分布函数f的计算公式如下：

式中，i表示粒子运动的方向，在D3Q19模型中，一共有19个方向，f_i ^eq表示平衡分布函数，ρ表示流体密度，e_i代表各个方向粒子运动速度的单位矢量，u表示动量，c＝Δx/Δt，Δx表示网格步长，Δt表示时间步长。计算过程如下：

1):由初始时刻t₀的ρ₀和u₀得到

2):通过迁移和碰撞得到新的粒子分布函数f_i ^eq；

3):计算新时刻t₁的ρ₁和u₁；

4):将ρ1和u1作为新的初始值进行下一时刻的计算。

LBM原始算法流程可查看图6。图7是本发明实施例中改进的LBM算法流程图：首先，将原算法的三个大循环拆分，按照腔内(体)、面、线(边)、点的顺序依次计算腔体内碰撞和迁移过程；其次，将计算网格点密度和质量部分单独取出，需要计算该网格点的质量时，将周围作用到该点的质量都算出来，然后再进行数据的迁移。由于一个向量处理单元一次可处理16个数据，因此，规定腔体中z轴方向(nz-2)为16的整数倍，其中，nz表示z轴的坐标长度。计算腔体内部时，因为不涉及边界，因此全是非固壁网格点的计算。然而，去除了边界部分，使得势能在x、y、z方向上连续，因此方便向量操作，可一次最多能计算16个网格点，大大提升了效率。由于之前规定z轴方向(nz-2)为16的整数倍，因此整个腔体可进行n次向量计算。

计算6个面的内部：同上，先去除4条边，先计算面内部的网格点，由于处于面上，其属于固壁边界，需要先进行壁面反射处理，即简单的数值交换，然后再计算网格点的数据碰撞和迁移。例如，当计算两个z平面时，由于z轴(nz-2)为16的整数倍，因此可以正好计算完成，计算其他平面时，本发明利用VLR向量长度寄存器相关函数。

计算12条线(即计算6个面的交线)时，通过VLR向量长度寄存器相关函数，设置好需要计算的向量单元长度，然后再进行交线上网格点的数据碰撞和迁移。

计算8个顶点(即计算各条线的交点)时，设置VLR长度为1，然后进行与线相同的操作，即通过VLR向量长度寄存器相关函数，设置好需要计算的向量单元长度，再进行各交点上网格点的数据碰撞和迁移。

至此，网格内所有点的势能均已计算完成。

同时，需要说明的是，在步骤S200中，设计面向基于向量的快速指数算法的具体过程如下：首先、针对求指数操作，将标准的泰勒展开式，改为n个连续的乘加操作；其次，设置多级系数字典，计算时通过查表取值；最后，利用指令延迟槽的循环展开。

以下举例说明：本发明针对exp操作设计将泰勒展开式，改为20级乘加的算法，除去第0级和第1级，共设置一个18个系数表，计算时通过查表取值值和利用指令延迟槽的循环展开，效率较标量算法提升了6*3*16倍。

exp运算的泰勒展开式如公式(5)，但在原有向量运算中，exp(x)运算实现起来效率低，所以设计了公式(6)，只需要利用乘加一种算术运算法则就能实现求解exp(x)，结构简单，有利于并行操作。

e^x＝((G₁x+G₂)x+G₃)x+G₄)x+.....+1)x+1 (6)

G_i＝1/(n-i+1)！ (7)

本发明根据上述改进的exp计算公式，设计将G_i取值i＝1,2.......19计算，并保存结果为一个字典，当执行exp运算时，只需要对应取值即可，所以设计后的公式只有乘加运算。公式中级数取20，能准确的计算结果，但是当x为负数时，泰勒公式需要取一个较大的n值才能保证exp计算的准确性，因此本发明计算x相反数的倒数，即exp(-x)＝1/exp(x)，从而规避了这一缺点。

针对上述exp运算的改进，本发明相应设计了利用指令延迟槽的循环展开算法，向量Load指令的执行时间是9拍，乘加运算执行6拍，向量store指令是4拍，本发明同时进行6个循环，如图7所示。exp函数从第1拍～第6拍执行Load指令读取x0～x5，第6拍开始执行公式(6)，利用x0～x5进行乘加运算，执行6*19次乘加运算之后，再执行6次store指令，将结果写入AM。参照图8可以得知(图8中，L代表Load指令，MA代表M1、M2、M3的3条并发的乘加指令，S代表store指令)，若顺序执行泰勒展开式的传统向量算法，6*3*16个exp(x)操作需要6*(9+19*6+4)＝762拍，而利用本循环展开算法，仅需要9+19*6+4+5＝132拍，效率提升约6倍。

此外，需要说明的是，前文使用改进的三维格子玻尔兹曼方法进行验证实例时，在计算6个面的内部过程中，描述的再计算网格点的数据碰撞和迁移具体表现为，计算每个网格的下一个时刻的互作用力时，通过缓存中间结果减少内存和外存之间的数据交换，复用格子的当前状态和下一时刻状态的存储空间，实现矩阵复用设计。在原算法中，为了保证每次进行消息传递时每个网格点上的数据不被覆盖，准备了两组数组进行存放。本发明根据LBM方法的计算特点，复用一套数组，提高缓存率，如图9所示，粗线内的网格为真实数据网格，虚线框为计算后的网格位置，在粗线的上、下、左、右分别扩展一行、两行、两列、一列网格，每次计算都只需要计算粗线内的网格点数据。以t时刻为例，将每一次计算完成后的新数值赋值给本网格点的左下角对应的网格点，比如当计算(i,j)网格点的数值时，更新后的数值存放在(i-1,j-1)处，这样之前的数据就不会被覆盖，与其相邻的几个网格点在计算时可以直接取值并得到正确的数值，为了避免计算碰撞时，覆盖尚未更新的格子的状态，本发明在计算x＝0和y＝0的格子的碰撞时，暂时把x＝nx-1和y＝ny-1的格子状态的缓存起来，等全部格子计算完成之后，再将缓存的数据更新对应的格子状态。对应的，在t+1时刻则存放在右上角对应的网格点处，同时，本发明设计在计算候分为奇数次和偶数次不同，必须按照此顺序进行：在t时刻是从左下角开始往右上角的顺序推进计算，在t+1时刻则相反，必须从右上角的顺序推进计算，拓展到三维也如此。

S300、基于单核的验证流程自适应到多核，实现快速并行计算；具体地，该步骤实现过程为：根据任务的大小和芯片核的个数，动态申请、自动分配DSP的数量，假设芯片核的个数为M，步骤如下：

Step1：根据芯片核的个数M，分别在x、y、z维度上分割为mx、my、mz个子腔体，满足mx×my×mz＝M；

Step2：计算a＝M/mx，b＝M/my，c＝M/mz，a,b,c分别为x、y、z维度上区间的长度；

Step3：针对x维度上长度为M的范围，分成了[0,a-1],[a,2a-1],......,[(mx-1)a,mx(a-1)]的mx个区间；每一个x区间对应的y维度上可划分[0,b-1],[b,2b-1],......,[(my-1)b,my(b-1)]的my个区间，每一个y区间对应z维度上可划分[0,c-1],[c,2c-1],......,[(mz-1)c,mz(c-1)]的mz个区间；需要说明的是，参照图10，以24个核，mx,my,mz分别取2,3,4为例，可以进一步理解此步骤；

Step4：根据Step3的区间划分将总任务分配到M个芯片核中进行计算。

上述算法实现了多核的自适应处理，能快速并行的算出结果。

总之，本发明针对多核向量程序的难点，通过预取、计算和写回三个状态的设计将DDR数据读写时间隐藏在DSP计算过程中，并设计基于表面切割的复杂标量程序的向量化方法，以及面向基于向量的快速指数算法，并使用改进的三维格子玻尔兹曼方法进行验证实例验证，修改内层循环，消除标量间的取模操作，实现效率提升，尤其，在针对指数运算执行效率不高的问题进行数学形式的转换，并给出多级系数字典用于计算时直接查表，最后，基于单核的验证流程实现自适应多核运算，能根据任务需求，自动分配每个DSP核的负载。相比现有技术，很好地解决了向量的并行运算的问题，缓解了DMA带宽低导致的时延长，实现了指数运算的高度向量并行计算。

上面的描述中阐述了很多具体细节以便于充分理解本发明，但是，本发明还可以采用其他不同于在此描述的其他方式来实施，因此，不能理解为对本发明保护范围的限制。

总之，本发明虽然列举了上述优选实施方式，但是应该说明，虽然本领域的技术人员可以进行各种变化和改型，除非这样的变化和改型偏离了本发明的范围，否则都应该包括在本发明的保护范围内。

Claims (6)

1.一种面向多核向量处理器的格子玻尔兹曼优化方法，其特征在于，包括以下步骤：

S100、通过设定的预取、计算和写回三个状态，针对芯片内核的向量访存部件进行数据传输时间隐藏的设计；

S200、设计基于表面切割的复杂标量程序的向量化方法，以及面向基于向量的快速指数算法，并使用改进的三维格子玻尔兹曼方法进行验证实例验证，实现效率提升；

S300、基于单核的验证流程自适应到多核，实现快速并行计算，该步骤的实现过程为：根据任务的大小和芯片核的个数，动态申请、自动分配DSP的数量，假设芯片核的个数为M，步骤如下：

Step1：根据芯片核的个数M，分别在x、y、z维度上分割为mx、my、mz个子腔体，满足mx×my×mz＝M；

Step2：计算a＝M/mx，b＝M/my，c＝M/mz，a,b,c分别为x、y、z维度上区间的长度；

Step3：针对x维度上长度为M的范围，分成了[0,a-1],[a,2a-1],......,[(mx-1)a,mx(a-1)]的mx个区间；每一个x区间对应的y维度上可划分[0,b-1],[b,2b-1],......,[(my-1)b,my(b-1)]的my个区间，每一个y区间对应z维度上可划分[0,c-1],[c,2c-1],......,[(mz-1)c,mz(c-1)]的mz个区间；

Step4：根据Step3的区间划分将总任务分配到M个芯片核中进行计算。

2.根据权利要求1所述的面向多核向量处理器的格子玻尔兹曼优化方法，其特征在于，所述步骤S100具体表现为：根据设定的预取、计算和写回三个状态依次轮换，将DDR数据读写时间隐藏在DSP计算过程，实现对芯片内核的向量访存部件的访存延迟隐藏。

3.根据权利要求2所述的面向多核向量处理器的格子玻尔兹曼优化方法，其特征在于，所述步骤S200中使用改进的三维格子玻尔兹曼方法进行验证实例具体按照体-面-边-点的维度分别进行向量计算，从而消除标量的取模操作。

4.根据权利要求3所述的面向多核向量处理器的格子玻尔兹曼优化方法，其特征在于，所述按照体-面-边-点的维度分别进行向量计算，从而消除标量的取模操作过程分解为如下步骤：

S201A、计算腔体内部，即计算非固壁网格点；

S202A、计算6个面的内部：先计算6个面内部的网格点，再计算网格点的数据碰撞和迁移；

S202A、计算6个面的交线：通过VLR向量长度寄存器相关函数，设置好需要计算的向量单元长度，再进行各交线上网格点的数据碰撞和迁移；

S204A、计算各条线的交点：设置VLR长度为1，然后通过VLR向量长度寄存器相关函数，设置好需要计算的向量单元长度，再进行各交点上网格点的数据碰撞和迁移。

5.根据权利要求4所述的面向多核向量处理器的格子玻尔兹曼优化方法，其特征在于，所述步骤S200中设计面向基于向量的快速指数算法的具体过程如下：

S201B、针对求指数操作，将标准的泰勒展开式，改为n个连续的乘加操作；

S202B、设置多级系数字典，计算时通过查表取值；

S203B、利用指令延迟槽的循环展开。

6.根据权利要求4所述的面向多核向量处理器的格子玻尔兹曼优化方法，其特征在于，所述步骤S202A中再计算网格点的数据碰撞和迁移具体表现为：

先计算每个网格的下一个时刻的互作用力时，通过缓存中间结果减少内存和外存之间的数据交换；

再复用网格的当前状态和下一时刻状态的存储空间，实现矩阵复用设计。

Images