CN102750150B - 基于x86架构的稠密矩阵乘法汇编代码自动生成方法 - Google Patents

Abstract

本发明公开了一种基于x86架构的稠密矩阵乘法汇编代码自动生成方法，属于计算机系统应用技术领域。本方法为：1)将稠密矩阵乘法中的源矩阵A、B分别划分为多个子块矩阵；其中，矩阵A划分的子块矩阵A_Mc*Kc大小为Mc*Kc，矩阵B划分的子块矩阵B_Kc*Nc大小为Kc*Nc；2)将每个子块矩阵A_Mc*Kc拷贝到连续的地址空间中，将每个子块矩阵B_Kc*Nc拷贝至连续的地址空间中；3)代码生成工具利用变化函数将输入的每一A_Mc*Kc*B_Kc*Nc子块矩阵乘法过程翻译为x86处理器支持的汇编代码。本发明的代码生成效率高，并可方便的移植到不同的x86处理器上。

Description

基于x86架构的稠密矩阵乘法汇编代码自动生成方法

技术领域

本发明涉及计算机系统应用技术领域。具体涉及一种面向x86体系架构实现的稠密矩阵乘法函数的汇编代码自动生成方法。

背景技术

代码自动生成方法已成为目前的一研究方向，许多应用领域都采用代码自动生成方法来降低开发周期，提高开发效率，比如专利ZL 200610152345.4，一种基于ASN类型对象的代码生成方法。

稠密矩阵乘法是基础线性子程序库(BLAS)中三级函数之一。它广泛应用于高性能计算领域，LINPACK测试结果严重依赖于稠密矩阵乘法函数的性能。BLAS三级中其余函数均可通过转换为稠密矩阵乘法来实现。许多CPU厂商自主研发了针对他们自己处理器系统架构特征的高性能基础线性代数子程序库，如Intel MKL和AMD ACML。

美国德克萨斯大学奥斯汀分校超级计算中心高性能计算组开发的基础线性代数子程序库(GotoBLAS)通过手工编写汇编代码实现稠密矩阵乘法函数。其性能与CPU厂商自主研发的基础线程代数子程序库(MKL，ACML)相当，但由于其(GotoBLAS)采用手工编写汇编代码来实现高效的矩阵乘法运算，这将导致稠密矩阵乘法函数的性能与处理器硬件结构特征联系紧密，比如高速缓存、物理寄存器、算术运算部件等。这些针对处理器硬件资源的优化方法往往具有平台特殊性，即针对某个特定处理器平台硬件资源种类和数量优化的稠密矩阵乘法汇编代码，往往能够在这个处理器平台上获得最优性能效果，但是当移植到另一台具有不同硬件资源的处理器平台上时，往往会出现程序性能下降的问题。GotoBLAS采取在不同的处理器上手工编写不同的稠密矩阵汇编代码来解决它在不同处理器平台间的可移植性问题。这是一种低效的代码开发方式，而且由于汇编代码的可读性比较差，代码的维护工作也需要领域专家来完成。

美国田纳西大学创新计算机实验室开发的自动调优线性代数软件(ATLAS)采用“软件的自动经验优化(AEOS)”方法，为给定的处理器“自动”选定最优的稠密矩阵实现方法。由于ATLAS生成的代码依赖底层编译器编译时实施的寄存器分配优化技术和指令调度优化技术，这两种优化技术往往对稠密矩阵乘法性能有巨大的影响。通用的底层编译器，如gcc，提供的优化技术往往需要适用于所有应用，因而不能像领域专家手工调优汇编代码一样，针对稠密矩阵乘法的算法特点充分利用物理寄存器资源和算术运算部件资源。因此，虽然ATLAS解决了稠密矩阵乘法函数开发周期长，可移植性差的问题，但是它生成的代码性能往往与手工优化编写的基础线性代数子程序库GotoBLAS以及CPU生产商自主研发的基础线性代数子程序库MKL，ACML存在5％～10％的性能差距。

发明内容

针对现有技术中手工编写稠密矩阵乘法函数的汇编代码开发效率低、可移植性差；和已有的自动调优矩阵乘法函数技术依赖底层编译器优化技术、性能低的缺点，本发明的目的在于提供一种基于x86架构的稠密矩阵乘法汇编代码生成方法，可针对不同x86处理器平台硬件资源的种类和数量自动生成稠密矩阵乘法函数的汇编代码。为降低底层编译器对生成代码性能的影响，本发明将传统编译器使用的优化技术：寄存器分配和指令调度，从编译器中“分离”出来，针对稠密矩阵乘法函数运算特点做出改进并给出实现，然后作用于本发明生成的汇编代码。

本发明解决的稠密矩阵乘法计算形式为：C＝A*B+C。在后续说明内容中，A和B为源操作矩阵，C为结果矩阵，其中矩阵A的大小表示为M*K，矩阵B的大小表示为K*N，矩阵C的大小表示为M*N。本发明采用的技术方案如下：

1)将矩阵A和B划分为较小的矩阵进行计算。首先将矩阵A按Mc*Kc的大小划分为多个子块，将矩阵B按Kc*Nc的大小划分为多个子块。然后逐次将每个子块矩阵A_Mc*Kc拷贝到连续的地址空间中，将每个子块矩阵B_Kc*Nc拷贝至连续的地址空间中。本发明主要解决将子块矩阵乘法的计算过程变为汇编代码实现问题，该子块矩阵乘法计算顺序采用通用的j-i-k迭代顺序，其中k为最内层迭代，j为最外层迭代。

2)利用已有代码生成工具(POET)，将1)中描述的子块矩阵乘法计算过程翻译为x86处理器支持的汇编代码。其中POET是一个可将源程序语言编写的代码根据翻译规则产生用户自定义程序语言类型代码的工具，源程序语言可以是C或者Fortran等。在本发明中，POET只作为一个代码翻译器工具，就像开发C程序使用的vc编辑器环境一样。“变换操作规则说明”在本发明中为开发人员根据通用优化技术实现的程序函数，本发明使用的变换操作包含：

a)对子块矩阵乘法计算过程实施循环分块、展开处理。变换操作过程描述如下，将循环迭代j和i的迭代步长分别变为N_r和M_r。使得在一次迭代k过程中，每次从第一个子块矩阵中读取M_r个元素，从第二个子块矩阵中读取N_r个元素，一个基本计算块将得到M_r*N_r个计算结果并与上一次循环迭代k计算结果进行累加，最后将整个k循环得到的结果累加到矩阵C中。对子块矩阵乘法计算过程实施循环展开处理，将循环k的迭代过程展开U_k次，可减少循环迭代开销。

b)对子块矩阵乘法计算过程中重复出现的数组访问模式替换为标量访问模式，减少冗余的从内存装载数据的操作次数，即将程序代码中所有出现的数组访问模式，如tmp＝A[i*Mr+k]*B[k*Nr+j]，替换为tmp＝la*lb，其中la＝A[i*Mr+k]，lb＝B[k*Nr+j]。

c)将经由a)-b)变换的子块矩阵乘法C语言代码逐句翻译为编译领域中常用的四元组中间表达式，表达式形式为(操作符，源操作数1，源操作数2，目的操作数)，其处理过程是使用C语言代码描述的。

d)对c)生成的四元组中间表达式通过合并，删除冗余等操作，生成向量化的四元组中间表达式，格式仍然为(操作符，源操作数1，源操作数2，目的操作数)。

e)将d)产生的向量化四元组表达式，根据x86指令手册逐一翻译为相应汇编代码，即将操作符更换为指令名，将操作数更换为真实汇编指令支持的描述格式。

3)上述处理过程即可在不同x86处理器平台上自动生成稠密矩阵汇编代码，缩短了传统手工编写稠密矩阵乘法汇编代码的开发周期，由于上述处理过程对处理器底层硬件资源涉及不多，因而可直接移植到不同的x86处理器上。与已有的自动调优矩阵乘法生成技术相比，本发明上述过程在目前公知的技术领域范围内是第一个能完全自动生成稠密矩阵乘法函数汇编代码的实例。为了使本发明生成的稠密矩阵汇编代码与手工编写的汇编代码有相当的性能效果，本发明对生成的汇编代码实施下述的优化操作：

a)结合子块矩阵乘法算法特点，改进传统编译器寄存器分配算法，将物理寄存器划分为不同组，并对乘法算法过程中出现的不同类型变量分配不同组的寄存器。

b)结合矩阵乘法算法特点，主要的程序基本块即为三层循环构成，因而应用传统编译器指令调度优化算法——“模调度”，对分配完寄存器的的汇编代码实施指令调度优化操作，该算法对程序中循环结构代码块的调度效果显著。

c)为提高子块矩阵乘法计算过程中数据访存效率，本发明在汇编代码中插入预取指令。由本发明提出的矩阵乘法算法可知，第二个子块矩阵按流水方式每一次迭代i、k的过程中访问不同的K_c*N_r的数据，然后与完整的第一个子块矩阵A_Mc*Kc相运算，得到M_c*N_r个计算结果。本发明在每次访问来自第i次循环读取的第二个子块矩阵的K_c*N_r的数据的代码中，提前对第i+1次循环要访问的第二个子块矩阵的K_c*N_r的数据-实施预取。并在每次两个子块矩阵的计算部分代码中，即循环迭代k的计算过程中，对第一个子块矩阵中的数据进行预取指令插入，实施预取。

相对于已有的自动代码生成工作，本发明在该步骤中采用的三种优化操作，均为独立于底层编译器(gcc)完整的程序实现，通过为用户提供不同的控制变量Sched，Pref(详见表1)，用来控制不同优化操作的实施。

4)上述步骤2)～3)即可生成高效的子块矩阵乘法汇编代码。由步骤1)中的分块算法描述可知，完整的稠密矩阵乘法计算过程，还包含将两个子块矩阵拷贝到连续的地址空间和中。为了进一步优化最终完整的矩阵乘法函数，本发明使用过程2)中同样过程将这两步拷贝操作翻译为汇编代码。由于拷贝操作对矩阵乘法整体性能影响较为微弱，在变为汇编代码后，不需要按步骤3)中同样过程实施优化。

与现有技术相比，本发明的积极效果为：

本发明提出了一种面向x86体系架构生成稠密矩阵乘法汇编代码的方法，并结合矩阵乘法算法的特点对生成的汇编代码采取了一系列优化措施。相对于传统手工优化的基础线性代数子程序库(GotoBLAS)，本发明提供了一种更为高效的稠密矩阵乘法函数开发方法，并可方便的移植到不同的x86处理器上。相对于已有的自动调优基础线性代数子程序库(ATLAS)，本发明提供了结合矩阵乘法算法特点的优化方案，能够更加高效的利用处理器硬件资源，降低了底层编译器对生成代码的影响。本发明定义了一系列可控变量，用以控制生成不同优化程度的矩阵乘法函数代码。该特性可以用来帮助分析不用优化技术对矩阵乘法函数性能影响的程度。同时增加本发明生成代码的灵活性和可扩展性。本发明在Intel Penryn，Nehalem和AMD Shanghai处理器上均能生成与CPU厂商自主研发的基础线性代数子程序库(MKL，ACML)性能相当的稠密矩阵乘法函数。除此之外，本发明在Intel最新发布的支持256位向量化指令的Sandy Bridge处理器上生成的稠密矩阵代码性能超过了IntelMKL0.8％。

附图说明

图1为矩阵乘法代码生成过程流程图。

图2一条计算程序语句的POET解析结果。

图3为四元组中间代码表达式生成流程图。

图4为向量化四元组表达式生成流程图。

图5为在Intel SandyBridge处理器上本发明生成的双精度稠密矩阵乘法函数与IntelMKL GotoBLAS和ATLAS的性能比较结果。

具体实施方式

下面结合附图对本发明的具体实施步骤进行描述，以便更好的理解本发明。

1、图1是本发明主要工作流程图。本发明基于手工调优实现的基础线性代数子程序库(GotoBLAS)中稠密矩阵乘法算法思想，首先将矩阵A划分为大小为Mc*Kc的子块矩阵，将矩阵B划分为大小为Kc*Nc的子块矩阵。由于在实际应用中参与计算的矩阵规模往往很大，只能存储于内存中，因而需要将矩阵分块，并将子块矩阵拷贝到连续的地址空间中来减少稠密矩阵乘法运算过程中访存失效和快表(TLB)失效的次数。由于所有子块矩阵乘法计算过程完全一样，因而本发明只生成子块矩阵乘法部分的汇编代码，矩阵分块和拷贝过程初始由手工编写C语言代码实现。

本发明中，矩阵A划分的子块矩阵规模为Mc*Kc，矩阵B划分的子块矩阵规模为Kc*Nc，其中Mc为矩阵A在行方向上的分块大小，Kc为矩阵A在列方向上的分块大小同时也是矩阵B在行方向上的分块大小，Nc为矩阵B在列方向上的分块大小。Mc，Nc，Kc的取值受限于处理器提供的高速缓存(cache)的容量，替换策略和访问延迟等因素。一般x86处理器都包含至少两级高速缓存结构，其中一级高速缓存访问速度快于二级高速缓存，但规模较小，往往不能放下整个子块矩阵计算需要访问的全部元素。矩阵分块的规模范围由如下两条原则限定，其中L₁代表一级高速缓存容量，L₂代表二级高速缓存容量，他们的大小可通过处理器手册获取：

①2*(K_r*N_c+M_c*N_r)+M_c*K_r<＝L₁

②M_c*K_c<＝L₂/2

本发明提供脚本测试程序以分块规模Mc，Nc，Kc作为传入参数，根据上述两条公式限定范围，自动搜索这三个变量在不同x86处理器上的最佳分块规模。

2、本发明的核心部分包括两个阶段。第一个阶段利用POET代码生成工具将子块矩阵乘法函数的C语言代码实现变换为汇编语言代码实现；第二个阶段对第一个阶段得到的汇编代码实现进行优化，得到与CPU厂商自主研发的基础线性代数库性能相当的稠密矩阵乘法代码。这两个阶段体现于图1中的“变换函数”和“优化技术”两个部分。下面结合图1分别对这两个阶段做进一步的描述。

第一个阶段：

子块矩阵乘法通用的计算过程主要包括j-i-k三层循环迭代，其中j为最外层循环索引，k为最内层循环索引。图1中POET代码生成器，接收子块矩阵乘法函数和变换函数。其中变换函数在本发明中为开发人员根据通用优化技术实现的程序函数。POET是一个比较成熟的代码翻译工具，本发明利用该工具对输入的子块矩阵乘法函数实施一系列变换操作函数，最终产生子块矩阵乘法汇编代码。该工具对于本发明的作用类似编程人员使用vc编辑器开发C语言程序的效用。发明过程并不对POET技术本身进行扩充与改进，发明过程中提到的变换函数是根据传统矩阵乘法优化实施步骤编写的符合POET语法规则的一系列程序函数，由本发明开发人员编写实现。需要发明人员对POET的使用有一定的学习基础。本发明通过提供给用户提供一组可调参数，便于用户对部分变换操作进行控制。本发明采用的变换操作包括下述五个步骤，分别是：

a)循环分块、展开操作。

首先对两个子块矩阵乘法计算过程中的循环j和循环i实施循环分块变换。本发明定义分块变量Nr代表循环j的迭代步长，定义分块变量Mr代表循环i的迭代步长。于是子块矩阵乘法最内层循环k的一次迭代计算过程变为：每次从第一个子块矩阵A_Mc*Kc中读取Mr个元素，从第二个子块矩阵B_Kc*Nc中读取Nr个元素，两两相乘后得到Nr*Mr个不同的结果。Mr和Nr的取值范围受限于处理器平台物理寄存器资源数量，一般该值可通过处理器手册获取。在不考虑处理器向量化特征的情况下，Mr和Nr应满足下述公式：Nr+Mr+Mr*Nr<R，其中R代表物理寄存器资源数量，有关向量化特征因素处理方案将在后续-内容中给出说明。

然后，为减少循环k的迭代开销，本发明定义循环展开因子Uk代表循环k的展开次数。由于循环展开次数过多会导致指令高速缓存溢出，会给子块矩阵乘法计算性能带来负影响，因而本发明设定Uk的上限为16，并将该变量作为一个可控变量，提供给用户，用户可以根据实际需要修改Uk的上界，也可以固定Uk的取值。

b)标量替换操作。

一般将操作tmp＝A[i*Mr+k]*B[k*Nr+j]称为数组访问模式，将变换操作la＝A[i*Mr+k]，lb＝B[k*Nr+j]，tmp＝la*lb称为标量替换变换。经过步骤1)操作后，每次循环k迭代计算过程中读取的第一个子块矩阵中的元素，如A[i*Mr+k]，都要在计算过程中与Nr个来自第二个子块矩阵的元素相乘一遍，即相同的子块矩阵元素要重复使用Nr次。同理，每个来自第二个子块矩阵的元素都要与Mr个来自第一个子矩阵的元素相乘一遍。在汇编代码技术中，一般先将数据从内存加载到寄存器，然后用寄存器中的数值进行计算会得到较高的执行效率。在编译原理中，由于不能保证数组访问模式是否访问同一个内存数据，因而每次都要将出现的数组访问模式数据先加载到一个寄存器中，再进行计算，因而在矩阵乘法计算过程中数组访问模式会产生大量冗余的数据加载操作，因为每个数组元素如A[i*Mr+k]都要重复使用Nr次。执行标量替换操作，可以由编程人员控制将数组元素赋值给一个标量，并用该标量替代在计算过程中该数组元素出现的所有位置。在编译原理中，这种处理方式可以保证在物理寄存器数量足够的情况下，每个数组元素只被加载到一个空闲的物理寄存器一次，该数组元素值会一直占用这个物理寄存器直到其生存周期结束。

本发明提供的变换步骤a)和b)是通用的优化C语言程序的优化技巧，为了便于某些用户只想得到优化的矩阵乘法的C语言函数代码，本发明向用户提供可控变量Asm，用来控制是否对当前得到的优化的C语言代码执行后续汇编代码生成变换操作过程，若用户设定Asm＝0，则从该步骤之后的全部变换操作都不再实施，本发明将生成一个优化的矩阵乘法C语言版本代码函数。

c)中间代码表达式变换操作。

传统编译原理技术通过词法分析，语法分析建立语法树，然后根据语法树生成中间代码表达式，常用的中间代码表达式分为四元式，三元式和间接三元式。本发明需要开发人员对POET工具的使用方式有一定的学习基础。POET对输入代码的语法解析结果相当于编译原理技术中构造的语法树。图2是POET对计算过程中一条语句的语法解析结果示例。其中，201～207为POET提供的语法解析结果句柄，使用人员可获得这些句柄(具体获取方式可学习POET使用语法得知)从而对语法树的结构进行重新组合。本发明利用POET对子块矩阵乘法代码语法解析结果将对应的C语言代码逐一翻译为四元组中间表达式，格式为(操作符，源操作数1，源操作数2，目的操作数)，即(opname,src1,src2,dest)。如果四元组中某个域没有值，则使用NA代替。图3为本发明提供的四元组中间代码表达式生成方法流程图。本发明将经由步骤1)和2)变换操作生成的子块矩阵乘法的C语言代码输入给POET代码生成器工具，POET依次读入C语言代码的每一条语句，判断当前读入语句是否为循环结构，即301表示的Next#结构。若当前读入的语句是循环结构，则将循环控制条件303表示的Loop#结构翻译为如304所示的一组四元组表达式，循环体部分由POET迭代解析判断。否则，表明当前读入的语句为一条表达式语句，如305所示。进一步判断当前表达式语句的右操作数是否为算数表达式，若是，则产生如307所示的计算四元组表达式，否则产生如308所示的赋值四元组表达式，有一个源操作数域为NA。本发明生成的子块矩阵乘法循环k部分四元组表达式形式如下所示：

1.(move,Kc,NA,k)

2.(label_k,NA,NA.NA)

3.(load,ba,offa,load1)

4.(load,bb,offb,load2)

5.(mul,load1,load2,tmp1)

6.(add,res1,tmp1,res1)

7.,,,,,,

8.(sub,k,1,k)

9.(jnez,k,NA,label_k)

其中，NA表示该部分没有值，offa表示从第一个子块矩阵中读取的元素在内存中的地址与拷贝空间首地址ba的偏移，offb表示从第二个子块矩阵中读取的元素在内存中的地址与拷贝空间首地址bb的偏移。

d)向量化变换操作。

向量化是目前众多x86处理器支持的一种高效优化技术通过使用单指令多操作数(SIMD)指令集，可以在一个指令执行过程中处理更多数据。图4为本发明采用的向量化四元组表达式生成过程流程图。为便于说明，设L代表向量指令所能操作的最大操作数位数，如L＝128，则表明一条SIMD指令可以操作128位的数据，这个数据可以是128位地址对齐的内存操作数也可以是128位向量寄存器中的数据。本发明为用户提供可控变量Vec_len，用来控制本发明最终生成128位的向量化汇编代码还是256位汇编代码。

设n代表矩阵乘法运算的数据类型所占的字节数(1字节＝8位)，如果矩阵乘法操作的数据类型为双精度浮点实数类型，则n＝8，一个浮点双精度操作数在内存中占8个字节。设E代表一个L长的向量寄存器中能装载的矩阵元素个数，则E＝L/(n*8)。在考虑向量化因素后，操作a)中循环分块Nr和Mr需要满足下述公式：③Mr/E+Nr+(Mr*Nr)/E<R。该公式为考虑向量化情况下生成汇编代码实施的寄存器分配条件。

本发明提供两种向量化代码技术处理方案对步骤c)得到的普通四元组表达式实施向量化变换操作，从而得到子块矩阵乘法的向量化四元组表达式。这两种方法在本发明中都只作用于第二个子矩阵元素的向量化操作上，本领域熟悉人员可知这些方法亦可作用于第一个子矩阵元素的向量化操作上。本发明为用户提供可控变量Shuf，取值范围为0/1，即用户可在变量设置文件里对该变量设定取值0或者1，用来控制生成不同的向量化处理方式的子块矩阵向量化代码，若Shuf＝0则使用第一种向量化方法，否则使用第二种向量化方法。其中变量设置文件由本发明提供，涵盖本说明中全部出现的“可控变量”，用来方便用户设置取值，其范围参见表1。若用户未设定固定取值，则本发明提供的搜索脚本程序可针对表1中每一种变量的每一种取值实施搜索。

Shuf＝0，使用装载并复制SIMD指令，如128位双精度浮点操作数指令为movddup，256位双精度浮点数操作指令为vbroadcast。这些向量化装载指令实现把从内存中读取的一个元素装载到L长的向量寄存器中，并复制E-1次到该寄存器的高位区域。因此，从第二个子矩阵中读取的Nr个元素仍然需要Nr次向量化装载操作。该操作对应图4中Shuf取0的分支操作。

Shuf＝1，使用洗牌SIMD指令，如128位双精度浮点操作数指令为pshufd，256位双精度浮点操作数指令为vpermilpd和vperm2f128。该指令实现将一个L长的向量寄存器中原本装载的E个矩阵元素的顺序进行交换或重新排列，使得E个元素在L中的每个位置上都出现一次，如(e1,e2)将被洗牌一次变为(e2,e1)。因此，从第二个子矩阵中读取的Nr个元素首先需要Nr/E次向量化装载操作，然后每个向量化装载结果都需要(E-1)次洗牌操作，装载和洗牌总操作次数满足Nr/E+Nr/E*(E-1)＝Nr次。该操作对应图4中shuf取1的分支操作。结合图4本发明的向量化方法实施步骤如下：

i.首先依次遍历c)中生成的每条四元组中间表达式，判断操作符类型是否为数据加载(load)操作；

ii.如果不是，当前四元组中间表达式将被保留，并将原始操作类型变换为对应的向量化表达方式，比如mul操作在128位SIMD中对应的向量化操作为mulpd；

iii.如果当前四元组中间表达式的操作类型为数据加载操作，则进一步判断该加载操作的元素是否为第一个子块矩阵中的元素：

1)如果是第一个子块矩阵中的元素并且该元素所在内存地址满足L长对齐，则保留当前四元组条目并将加载操作load替换为向量化数据装载操作；

2)否则删除当前四元组中间表达式条目并递归删除所有使用该装载操作结果的计算四元组表达式条目；

iv.如果是第二个子块矩阵中的元素，则判断当前Shuf的取值是0——使用装载并复制向量化操作还是1——使用装载洗牌向量化操作：

1.如果是装载并复制向量化操作则替换第二个子矩阵原始数据加载操作为向量化数据装载并复制操作，对应图4中Shuf＝0分支操作；

2.如果是装载洗牌向量化操作则先将当前数据装载条目更换为向量化数据装载操作，然后删除第二个子矩阵后续E-1次矩阵元素装载操作，最后插入E-1次对当前向量化装载结果的洗牌操作，对应图4中Shuf＝1的分支操作。

e)汇编指令映射变换操作。

最后一步是将向量化的四元组表达式翻译为x86处理器支持的汇编指令。由于矩阵乘法操作类型非常有限，主要包括数据装载操作，乘法操作，加法操作，洗牌操作和写回操作，因此可以根据x86处理器指令手册规定不同四元组表达式映射结果。比如当前要映射的向量化四元组表达式为(mulpd,load1,load2,tmp1)，完成tmp1＝load1*load2计算。由于x86处理器支持的128位SIMD乘法指令只有三个域，本发明将(mulpd,load1,load2,tmp1)翻译为(movapd,load2,tmp2),(mulpd,load1,load2)两条指令，其中第一条指令用于保存load2的原始数值，第二条指令计算load1*load2并将计算结果存于load2中。然后，遍历后续四元组表达式，将原本读取tmp1的四元组表达式更换为读取load2，同时将原本读取load2的四元组表达式更换为读取tmp2。

下述代码片段为本发明第一阶段最终生成的x86处理器支持的汇编代码示例，在该实例中，矩阵乘法元素类型为浮点双精度类型，生成汇编代码为64位模式：

1.movq Kc,k

2.label_k:

3.movapd offa*8(ba),load1

4.movapd offb*8(bb),load2

5.movapd load2,tmp2

6.mulpd load1,load2

7.addpd load2,res1

8.,,,,,,

9.subq 1,k

10.jne k,label_k

第二个阶段

3、在第一个阶段中，本发明利用POET代码生成器，对子块矩阵乘法函数实施五步变换操作，将其C语言实现代码变换为汇编代码。至此，本发明即可在不同x86处理器平台上自动完成稠密矩阵汇编代码生成工作，相对于传统手工编写稠密矩阵乘法汇编开发过程，本发明可大大缩短代码的开发周期。此外，由于上述处理过程对处理器底层硬件资源涉及不多，因而可直接移植到不同的x86处理器上。为了使本发明生成的矩阵乘法汇编代码能够获得与CPU厂商自主研发的基础线性代数库MKL和ACML相当的性能效果，本发明在第二个阶段对之前第一个阶段生成的子块矩阵乘法汇编代码实施三步优化方案，这三个优化技术方案均为独立于底层编译器(gcc)由本发明完整实现的程序函数，具体技术方案实施步骤描述如下：

a)寄存器分配优化操作。

已有的稠密矩阵乘法自动调优技术，如ATLAS，往往依赖底层编译器(gcc等)对自动生成的代码实施寄存器分配优化技术。由于gcc是一种通用的编译器，它不能针对矩阵乘法算法特点对计算过程中出现的变量合理分配使用处理器提供的寄存器资源。因此相对于手工优化的矩阵乘法汇编代码而言，它的优化结果往往是不充分的。为了充分利用寄存器资源并减少由寄存器分配产生的数据假依赖关系，本发明提出一种结合矩阵乘法算法特点的寄存器分配方法。该方法将向量寄存器划分为四组：

i.第一组用来装载每次循环迭代读取的Mr个第一个子矩阵中的元素，初始值为Mr/E个向量寄存器；

ii.第二组用来装载每次循环迭代读取的Nr个第二个子矩阵中的元素，初始值为Nr/E个向量寄存器；

iii.第三组用来保存计算得到的Mr*Nr个结果，初始值为(Mr*Nr)/E个向量寄存器；

iv.第四组用来保存计算过程中设定要保存的临时结果数值，初始值为R-Mr/E-Nr-Mr*Nr/E个向量寄存器，其中R为处理器含有的总物理向量寄存器数量，这些临时结果由算法确定。

这四组寄存器都采用队列的方式实施申请和释放操作。每一种变量类型从其对应的寄存器队列头部申请一个空闲的寄存器，使用完毕后将释放的寄存器排于队尾。这种方式可以保证每个寄存器的使用情况相对平均，不会出现一直空闲或一直被占用的寄存器。寄存器分配策略在子块矩阵乘法的一开始就先将(Mr*Nr)/E个寄存器分配给结果数据，这(Mr*Nr)/E个寄存器会在一次Nr-Mr-Kc迭代结束后将保存的结果数据累加到结果矩阵C中再释放。其余三组寄存器按之前分配的初始值实施分配过程，如果在给子块矩阵计算过程中出现的变量分配寄存器时有一个寄存器组出现寄存器不够现象(即溢出)，本发明将从另外两组寄存器中找出寄存器数量最多的一组，取出一个寄存器给当前导致溢出的寄存器组。然后重新实施寄存器分配过程，直至整个子块矩阵计算过程中出现的变量都被合理分配完寄存器为止。

b)由于主流x86处理器普遍支持超标量指令发射技术，即在同一个时钟周期内可以同时执行多条互相独立的指令操作，因此一个精心排布的指令序列往往能获得更好的性能效果。指令调度，是将汇编代码中没有依赖关系的指令放在一起，以便在同一个指令周期被处理器读取执行。已有的自动调优稠密矩阵代码生成技术很少有对生成的代码实施指令调度优化方案，如果依赖底层编译器对生成的代码实施指令调度优化操作，由于编译器在寄存器分配阶段的优化不够充分，将直接影响代码潜在并行性效果。

子块矩阵乘法运算过程中的循环k程序块是整个代码的热点区域，如果该代码片段的指令排布合理，能充分利用处理器的计算部件，则整个子块矩阵乘法函数往往可以取得高效的性能效果。因此，本发明采用传统的模调度(Modulo Scheduling)算法对子块矩阵乘法最内层循环k程序块实施指令重排序。该调度算法为一个成熟的技术，相关实施操作方法可在各大论文数据库中检索得到，而一般程序员皆可通过阅读相关论文资料，对该算法提出的实施方案给出具体的程序代码实现。

本发明对该调度算法给出了一个C语言程序实现，并为用户提供可控变量Sched，控制是否将该优化操作技术作用于本发明生成的子块矩阵乘法汇编代码上，若Sched＝0，则表明不对生成的汇编代码实施指令优化操作。该功能可以控制生成一个不带指令调度的子块矩阵乘法汇编代码，可以方便领域熟悉人员在此基础上自行实现手工代码调度优化方法，亦可以作为扩展接口，实现其他指令调度方案，将其作为新的调度组件加入本发明中。

该指令调度算法本身并非本发明提出的新技术，相对于已有的自动调优矩阵乘法代码生成技术，本发明并非依赖底层编译器对生成的代码实施指令调度优化操作，而是提供了一种更适合矩阵乘法特点的循环调度优化技术，即将“模调度”指令调度算法通过编程实现，作为一个优化模块融合到本发明代码生成机制中，由用户控制对特定的代码热点，循环迭代k，实施指令调度优化，是本发明所提出的新的技术解决方案。

c)预取指令可以把数据在实际使用前从内存移动到高速缓存(cache)中，从而将访存的失效开销与数据的计算时间相重叠以提高程序的执行效率。众所周知，数据在高速缓存中被处理器执行的速度远远快于数据在内存中的情况。当今x86处理器普遍支持多层高速缓存结构，一级高速缓存访问速度要高于二级甚至三级高速缓存。由于实际应用中矩阵A，B，C的规模往往比较大，一般都会存储在内存中，为了提高数据重用性和空间局部性，本发明在实施步骤1中首先将矩阵A，B，C换分为小块，然后拷贝至两个连续的地址空间中。由于一级高速缓存规模比较小，通常不能放下完整的两个子块矩阵规模，因而划分的两个子块矩阵一般存储于二级高速缓存中，在矩阵乘法计算过程中通过读取(load)数据操作，一次从二级或三级高速缓存中加载一个高速缓存数据块(cache line)到一级高速缓存中。

矩阵乘法运算过程比较规则，子块矩阵的循环迭代顺序为j-i-k，经过循环分块和展开操作后，j每迭代一次，第二个子块矩阵B将读取不同的Kc*Nr个数据，然后与整个Mc*Kc的第一个子块矩阵相乘。因此，本发明采用在每次i迭代前对下一次要读取的Nr*Kc的第二个子块矩阵元素实施预取操作，即在此位置插入预取指令。然后i每迭代一次，第一个子块矩阵A将读取不同的Mr*Kc个数据，于是在每次Mr*Kc的迭代计算过程中对下一次要读取的Mr*Kc个数据实施预取操作，即在k开始迭代前插入对第二个子块矩阵A的预取指令。结果矩阵C在一次j-i迭代中只需要被访问Mr*Nr个数据元素，访问数据量较少。另一方面，由于未对结果矩阵C实施分块拷贝处理，结果矩阵C的数据元素仍然存储在内存中，数据局部性很差，本发明在每次i迭代开始前对当前计算结果要写回的结果矩阵C中的元素实施预取操作，这样，在k迭代结束后，要访问的结果矩阵C中的元素应该已经提前加载到一级高速缓存中。

x86处理器支持4种数据预取指令，prefetcht0，prefetcht1，prefetcht2，prefetchnta，用以将内存中的数据提前移动到不同高速缓存层次中，具体指令含义可在x86指令手册中查到。在本发明中对两个子块矩阵以及结果矩阵C都有数据预取操作，但将这三类数据预取到哪一级高速缓存中，在不同的处理器上往往会选取不同的预取指令实现。本发明为了能够在各个x86处理器平台上均获得最佳预取效果，为用户提供了可调变量Pref_a，Pref_b，Pref_c来替代实际汇编代码中预取指令的名字。比如实际子块矩阵乘法汇编中按上述方法和位置插入对第一个子块矩阵的预取指令为Pref_a然后在汇编代码的开头对Pref_a实施宏替换，即可将Pref_a映射为上述提到的四种不同的预取指令。如果用户为这些变量语句设置不同的宏值，则可实现不同的预取指令预取数据的效果。本发明在预取优化提出的新的处理技术是使用宏定义来表示预取指令，然后由用户修改宏定义来尝试不同的预取指令效果，从而可得到最佳的预取指令组合结果。

在AMD处理器上一个实施预取指令的经验是使用prefetchw指令。该指令只在AMD处理平台上获得支持，该指令对将会写回缓存的数据有显著的预取效果。在矩阵乘法中，结果矩阵的数据访问符合这种情况，因此在AMD处理器上，本发明固定Pref_c的取值为prefetchw。为了给某些领域熟悉人员提供手动插入汇编指令的需要，本发明为用户提供了可控变量Pref，控制该优化操作是否实施。若Pref＝0，则不对子矩阵乘法汇编代码插入数据预取指令。否则，本发明会调用预取插入程序，对生成汇编代码按上述方法插入对两个子块矩阵以及结果矩阵的预取指令。

4、在得到子块矩阵乘法函数的汇编代码后，本发明为降低步骤1中矩阵乘法分块拷贝操作的开销，进一步提高整体矩阵乘法函数的性能，利用第一阶段的发明过程将拷贝操作也自动转换为汇编代码操作。由于拷贝操作中仅包含数据的装载和写回操作，并且写回的地址为连续对齐的空间地址，因此拷贝操作是矩阵乘法操作的一个子集，可以直接由本发明提供的代码生成优化技术生成汇编代码。

表1是本发明过程中定义的全部变量，分为与处理器硬件资源相关变量和代码生成过程控制变量。这些变量可由使用本发明的用户指定具体定值，也可在未设定定值的情况下，根据表1中给出的取值范围，由本发明提供的脚本测试程序，自动测试所有可能取值结果找到适用于给定x86处理器平台最优的矩阵乘法汇编代码。本发明提供一个脚本程序来完成上述所说的自动测试过程，该脚本程序可以根据表1中与处理器硬件资源相关变量可能的不同取值生成不同的子块矩阵乘法汇编代码以及拷贝操作汇编代码，并使用2048*2048的矩阵规模对完整的矩阵乘法函数进行性能测试。本发明将测试结果以及生成相应汇编代码的变量取值一同记录到一个名为FINAL文件中。在尝试完所有变量取值之后，本发明检索FINAL文档，找出性能最好的结果以及生成该结果的变量的取值，然后根据这些取值重新生成最终完整的矩阵乘法函数。也可以由有经验的开发人员自行设定经验最优值，当需要自动测试找到最优代码参数时，表1提供了传统优化经验的参数取值范围，其中Mc，Nc和Kc参考公式①,②获得取值范围，Nr和Mr参考公式③获得取值范围。

表1、变量列表

注：表1中[a:b:c]格式表示下界取值为a，上界取值为b，取值间隔为c。

显然，本领域的技术人员应该明白，上述的本发明的各个步骤可以移植到非x86体系架构处理器上。只需要要更改步骤2中变换处理e)将向量化的四元组表达式映射到其它处理器架构上支持的指令集即可。这样，本发明不局限于x86处理器本身。凡在本发明的思想和原则内，作出的任何修改，等同变换，改进等均应该包含在本发明的保护范围之内。

Claims (12)

1.一种基于x86架构的稠密矩阵乘法汇编代码自动生成方法，其步骤为：

1)将稠密矩阵乘法中的源矩阵A、B分别划分为多个子块矩阵；其中，矩阵A划分的子块矩阵A_Mc*Kc大小为Mc*Kc，矩阵B划分的子块矩阵B_Kc*Nc大小为Kc*Nc；Mc为矩阵A在行方向上的分块大小，Nc为矩阵B在列方向上的分块大小，Kc为矩阵A在列方向上的分块大小，矩阵A列值与矩阵B行值相等；

2)将每个子块矩阵A_Mc*Kc拷贝到连续的地址空间中，将每个子块矩阵B_Kc*Nc拷贝至连续的地址空间中；

3)代码生成工具利用变化函数将输入的每一A_Mc*Kc*B_Kc*Nc子块矩阵乘法过程翻译为x86处理器支持的汇编代码。

2.如权利要求1所述的方法，其特征在于代码生成工具生成所述汇编代码的方法为：

a)对输入的A_Mc*Kc*B_Kc*Nc子块矩阵乘法进行循环分块、展开操作：在每一次迭代k过程中，每次从第一个子块矩阵A_Mc*Kc中读取Mr个元素，从第二个子块矩阵B_Kc*Nc中读取Nr个元素，两两相乘后得到Nr*Mr个结果并将其累加到一结果矩阵C中；其中，Nr代表所述稠密矩阵乘法的j-i-k三层循环迭代计算过程中循环j的迭代步长，Mr代表所述稠密矩阵乘法的j-i-k三层循环迭代计算过程中循环i的迭代步长；

b)中间代码表达式变换操作：将步骤a)的子块矩阵乘法源程序语言代码逐句翻译为四元组中间表达式，表达式形式为(操作符，源操作数1，源操作数2，目的操作数)；

c)向量化变换操作：对所生成的四元组中间表达式进行合并、删除冗余处理，生成向量化的四元组中间表达式；

d)汇编指令映射变换操作：将向量化四元组表达式逐一翻译为相应汇编代码。

3.如权利要求2所述的方法，其特征在于设定一循环展开因子Uk，用于控制循环k的展开次数。

4.如权利要求2所述的方法，其特征在于对子块矩阵乘法计算过程中，将重复出现的数组访问模式替换为标量访问模式。

5.如权利要求2所述的方法，其特征在于将步骤a)的子块矩阵乘法源程序语言代码逐句翻译为四元组中间表达式的方法为：代码生成工具依次读入源程序语言代码的每一条语句，判断当前读入语句是否为循环结构，如果是循环结构，则将循环控制条件结构翻译为一组四元组表达式，否则当前读入的语句为一条表达式语句；进一步判断当前表达式语句的右操作数是否为算数表达式，若是，则产生一计算四元组表达式，否则产生一赋值四元组表达式。

6.如权利要求2所述的方法，其特征在于生成向量化的四元组中间表达式的方法为：

1)依次遍历每条四元组中间表达式，如果其操作符类型不是数据加载操作，则保留当前四元组中间表达式，并将原始操作类型变换为对应的向量化表达方式；

2)如果当前四元组中间表达式的操作类型为数据加载操作：

21)如果加载操作的元素是第一个子块矩阵中的元素，且该元素所在内存地址满足L长对齐，则保留当前四元组条目并将加载操作替换为向量化数据装载操作；如果该元素所在内存地址不满足L长对齐，则删除当前四元组中间表达式条目并递归删除所有使用该装载操作结果的计算四元组表达式条目；

22)如果加载操作的元素是第二个子块矩阵中的元素，且该元素使用装载并复制向量化操作，则替换第二个子矩阵原始数据加载操作为向量化数据装载并复制操作；如果该元素使用装载洗牌向量化操作，则先将当前数据装载条目更换为向量化数据装载操作，然后删除第二个子矩阵后续E-1次矩阵元素装载操作，最后插入E-1次对当前向量化装载结果的洗牌操作：

其中，L为装载所读取元素的向量寄存器的长度，E代表一个L长的向量寄存器中能装载的矩阵元素个数，Nr和Mr满足关系式：Mr/E+Nr+(Mr*Nr)/E<R，R代表物理寄存器资源数量。

7.如权利要求2所述的方法，其特征在于Nr和Mr的取值满足关系式：Nr+Mr+Mr*Nr<R；其中R代表物理寄存器资源数量。

8.如权利要求2所述的方法，其特征在于所述源程序语言为C语言或Fortran语言。

9.如权利要求2所述的方法，其特征在于对所述汇编代码进行调优，其方法为：将向量寄存器划分为四组，第一组用来装载每次循环迭代读取的Mr个第一个子块矩阵中的元素，初始值为Mr/E个向量寄存器；第二组用来装载每次循环迭代读取的Nr个第二个子矩阵中的元素，初始值为Nr/E个向量寄存器；第三组用来保存计算得到的Mr*Nr个结果，初始值为(Mr*Nr)/E个向量寄存器；第四组用来保存设定的计算过程中需要临时保存的数值，初始值为R-Mr/E-Nr-Mr*Nr/E个向量寄存器；其中，L为向量寄存器的长度，E代表一个L长的向量寄存器中能装载的矩阵元素个数，R为总物理向量寄存器数量。

10.如权利要求9所述的方法，其特征在于四组所述向量寄存器均采用队列的方式实施申请和释放操作；如果第一、二、四组向量寄存器出现寄存器溢出时，从另外两组寄存器中找出寄存器数量最多的一组，取出一个寄存器给当前导致溢出的寄存器组，然后重新实施寄存器分配，直至整个子块矩阵计算过程中出现的变量都被合理分配完寄存器为止。

11.如权利要求1所述的方法，其特征在于对所述汇编代码进行调优，其方法为：采用模调度算法对子块矩阵乘法最内层循环k程序块实施指令重排序。

12.如权利要求2所述的方法，其特征在于对所述汇编代码进行调优，其方法为：步骤a)中，当每次访问来自第i次循环读取的第二个子块矩阵的K_c*N_r的数据的代码时，提前对第i+1次循环要访问的第二个子块矩阵的K_c*N_r的数据实施预取当每次访问来自第二个子块矩阵的K_c*N_r的数据的代码时，对下一次要访问的K_c*N_r的数据插入预取代码，并在每次两个子块矩阵的计算部分代码中对第一个子块矩阵中的数据实施预取。