CN102722472B - 一种复数矩阵的优化方法 - Google Patents

Abstract

本发明公开了一种复数矩阵的优化方法，特征是先计算在龙芯体系结构下特定的展开粒度，对复数矩阵进行四乘四循环展开，将矩阵分块大小nb选取为一个小于52并使得24乘以矩阵分块大小nb的平方小于龙芯处理器一级数据缓存64千字节的最大的值，得到在龙芯上最优的理想矩阵分块大小；并利用数据存储的连续性和局部性对矩阵乘法中的矩阵进行合理划分与合并，减少龙芯的一级数据缓存访问次数；对复数矩阵运算中的普通复数乘加运算利用复数经典算法中两个复数之间的乘法算法，从而减小了运行规模，使得复数矩阵乘法运算在龙芯3A上的计算性能提高了大约50％，可使BLAS库在龙芯3A上的运算速度达到原来的2.5倍以上。

Description

一种复数矩阵的优化方法

技术领域

本发明属于计算机程序优化方法技术领域，具体涉及中函数优化设计复数矩阵函数的优化方法。

背景技术

随着大规模数学计算在科学、工程、金融等领域的广泛应用，常用数学库基础线性代数程序集(BLAS)已成为各种系统软件必需的矩阵运算工具。另外，对高性能计算进行评价的重要测试线性系统软件包(LINPACK)正是基于基础线性代数程序(BLAS)库进行测试，因此，在不同平台上对基础线性代数程序(BLAS)库的优化是非常重要的。龙芯3A是中国科学院计算技术研究所研制的首款国产四核高性能处理器，采用65nm工艺制造，最高工作主频为1GHz，该芯片内集成了四个64位的四发射超标量GS464高性能处理器核，每个处理器核包含64KB的一级数据缓存(cache)和64KB指令缓存，芯片内集成四核共享的4MB二级Cache，采用128位接口的交叉开关网络(AXI)。其浮点性能的好坏与BLAS库的性能密切相关。而BLAS库中约一半的函数是复数浮点运算函数，因此，对复数矩阵运算进行优化十分必要。以往的类似优化都是针对通用矩阵或特定的其他处理器平台进行优化。基础线性代数程序(BLAS)库本身并没有为复数设计单独的数据结构来完成库函数的运算，导致针对复数的运算规模较大，整体效率较低，从而运行开销较大。

发明内容

本发明的目的是提出一种复数矩阵的优化方法，以克服已有通用的整块进行计算和固定分块大小的计算方法在龙芯上运行时整体效率较低而导致运行开销大的问题。

本发明在龙芯上的复数矩阵的优化方法，其特征在于操作步骤如下：

第一步：利用展开粒度公式K_AK_B/(K_A+K_B)＝f/g计算最优取数比，其中K_A和K_B分别为矩阵A和矩阵B的展开粒度，f/g为最优运算取数比，令两个矩阵展开粒度的乘积除以两个矩阵的展开粒度之和所得的值与最优取数比相同，在龙芯上调整两个矩阵的展开粒度K_A和K_B，使得在循环结构中保证矩阵A和矩阵B的运算取数比f/g为2：1，则所得的展开粒度即为在运算中一次循环需要展开的运算次数；

第二步：在矩阵运算时，对于矩阵内数据大于500的大规模矩阵计算，对每一分块进行运算时在循环中采用4×4的展开，每个展开的运算对应龙芯的一个发射槽，通过发射槽发射运算指令给龙芯运算器进行乘加运算；根据矩阵分块公式24×(nb)²≤NK，其中NK为龙芯的一级数据缓存(cache)的大小，将矩阵分块大小nb选取为一个小于52并使得24乘以矩阵分块大小nb的平方(nb)²小于龙芯处理器一级数据缓存64千字节(kb)的最大的值，得到在龙芯上最优的理想矩阵分块大小；

第三步：对于在复数矩阵运算中出现的普通复数乘加运算，利用复数经典算法中两个复数之间的乘法算法，设第一复数为A＝x+yi，第二复数为B＝m+ni，第三复数为C＝p+qi,计算第一复数A乘以第二复数B后再加第三复数C的值，将此值赋给第三复数作为第三复数的新值，算法流程如下：计算第一复数A的实部x乘以第二复数B的实部m的值为第一实数e(e←x×m)，计算第一复数A的虚部y乘以第二复数B的虚部n的值为第二实数f(f←y×n)，计算第三复数C的实部p加第一实数e后减去第二实数f的值，并将此值赋给第三复数C的实部p(p←p+e-f)，计算第一复数A的实部x加第一复数A的虚部y的值为第三实数k(k←x+y),计算第二复数B的实部m加第二复数B的虚部n的值为第四实数t(t←m+n)，计算第三复数C的虚部q加第三实数k乘以第四实数t的积再减去第一实数e和第二实数f的值，并将此值赋给第三实数k的虚部q(q←q+k×t-e-f)，则得到第三复数C的新值，通过龙芯发射槽传给运算器参与分块矩阵的运算。

由于现有技术采取的是直接将矩阵全部带入进行计算，或利用通用的固定分块方法进行计算，若在龙芯上用现有技术进行优化运算，就会产生问题规模较大、整体效率较低而导致运行开销大的问题。而本发明的复数矩阵的优化方法，是针对龙芯3A的四发射体系结构和特定的一级数据缓存大小进行的优化，通过计算特定的展开粒度，对复数矩阵进行四乘四循环展开，特定大小分块，并对复数矩阵运算中的普通复数乘加运算利用复数经典算法减小了运行规模，进而使得复数矩阵乘法运算在龙芯3A上的计算性能提高了大约50％，而由于在BLAS库中约一半的函数是复数矩阵浮点运算函数，因此采用本发明复数矩阵的优化方法可使BLAS库在龙芯3A上的运算速度达到原来的2.5倍以上。提高了整体运算效率。

附图说明

图1是在龙芯上改变使用不同分块大小后常用数学库基础线性代数程序集(BLAS)中一个名为zGEMM的广义矩阵相乘的函数的性能测试比较的矩阵分块性能图。

图2为采用本发明复数矩阵的优化方法在龙芯处理器上的优化效果仿真图。

具体实施方式

实施例1：

本发明在龙芯上的复数矩阵的优化方法实施例的具体操作步骤如下：

1、循环展开和指令调度

常用数学库基础线性代数程序集(BLAS)中的三级BLAS的循环展开一般指的是矩阵分块后，块状矩阵做乘法运算时的展开粒度，侧重考虑循环展开和指令调度之间的关系。由于龙芯四发射的特性：每个周期发射一条取数指令，两条浮点指令，一条定点指令，为了尽可能缩短程序运行的周期，就要在循环结构中保证运算取数比为2：1，设两个矩阵第一矩阵A和第二矩阵B的展开粒度分别为KA和KB，指令调度时就需要满足KA KB/(KA+KB)＝2，且KA和KB为小于矩阵分块大小的正整数。在该范围内很容易穷举出只有(4，4)和(3，6)组合，在展开度相差不大的情况下，可读性更强、访存一致性更好的四乘四展开就是最优的选择。

2、矩阵分块

大规模计算时，为了有效减少龙芯处理器的一级数据缓存Cache的失效，会使用矩阵分块技术，对矩阵进行合理分块，分块大小对浮点运算性能的影响很大，实验表明，每次参与运算的矩阵块大小保证在龙芯一级数据缓存Cache以内为最优。设矩阵分块大小为nb，zGEMM在块内实现的是两个矩阵分块大小nb×nb的乘法，另外由于每次还要保存一个中间累加的结果，故参与运算的数据大小为3×8×nb×nb＝24×nb×nb个字节bytes，龙芯处理器的一级数据缓存cache大小为64KB，为了保证不会产生容量失效，需要保证24×nb×nb不大于64K，同时由于4路组相联映射本身会有自冲突失效的可能性，nb的选取应该比计算得到的值要小。此外若矩阵分块大小nb太小的话，浮点性能峰值会很低，因为过多的分块会造成访存次数的增多和分支转移预测成功率的下降，不能充分发挥多级浮点指令流水线的性能，龙芯处理器上zGEMM优化用的是四乘四循环展开，当24nb×nb≤64K且nb为4倍数时nb＝52，所以实际应用中矩阵分块大小nb选择的是一个小于52且足够大的值。

图1是在龙芯上改变使用不同分块大小后常用数学库基础线性代数程序集(BLAS)中一个广义矩阵相乘的函数zGEMM的性能测试比较的矩阵分块性能图。由图1可见：最上面的曲线48表示矩阵分块大小nb＝48时的运算速度，此时的性能最高，然后依次次是矩阵分块大小为56和40的表示运算速度的曲线，矩阵分块大小为36时的运算速度的曲线在最下面，表示其性能最差。由图1中的曲线说明了矩阵分块大小为48是zGEMM一个比较好的分块大小，这和理论分析的结果比较吻合。在优化方法一定的情况下，它在访存次数过多和容量失效之间寻找到了一个有效的平衡点。利用数据存储的连续性和局部性对矩阵乘法中的矩阵进行合理划分与合并，减少龙芯一级数据缓存cache访问次数，降低了计算的规模，从而大幅提高了BLAS库复数运算在龙芯3A上的运算效率。

3、运算访存比

三级BLAS优化过程中也需要充分考虑运算访存比的影响，在使用循环展开、指令调度和矩阵分块后，矩阵分块大小nb×nb的基本块实现时，每一轮内循环完成对结果矩阵C的4×4矩阵块的更新，4×4小矩阵块更新过程中涉及到的访存量分别为：Ma＝4×nb，Mb＝4×nb，Mc＝4×4×2，Mabc＝8×nb+32；涉及到的运算量为：使用乘加指令后I＝16×nb次乘加，龙芯处理器每个周期cycle可以同时发射一条访存指令，两条浮点指令，因此不考虑其它因素时，每一轮更新的4×4小矩阵块的访存的周期数为Mabc，浮点操作的周期数为I/2，理论上的效率可以达到：Mabc/(I/2)＝1+4/nb，由此可见分块过大的话也会降低运算访存比，分块大小要使得运算访存比适中。利用数据的局部性对复数运算进行上层合并以减少访龙芯一级数据缓存cache次数，进而大大提高了复数乘法运算在龙芯3A上的计算性能。

4、复数乘法优化

由于矩阵分块的存在，zGEMM完成规模为N方阵运算时，要进行个分块矩阵乘法和个分块矩阵加法的操作，也就是说块内实现的是规模为矩阵分块大小nb的4乘4加的矩阵运算。

三级BLAS函数完成矩阵和矩阵的运算，数据的运算量是O(n3)量级。在处理大规模矩阵计算时CPU时钟周期更多消耗在浮点计算中，故属于浮点计算密集型函数，它性能提升的空间很大，如何更高效地开发处理器浮点运算潜能，是该类函数的优化目标。三级BLAS复数运算函数由于运算量非常大，一般是实数运算函数的4倍，为了避免更多的访存失效，BLAS会对输入的复数矩阵进行实部和虚部的分离，然后再调用同类型的实数运算函数来完成分离后矩阵间的运算，最后再将计算结果整合在一起。

复数分离后的实部和虚部矩阵是连续存储的，利用数据的局部性可以对复数运算进行上层合并从而减少访cache次数，这样运算取数比从2：1可以增加到2.4：1，理论上有一定性能的提升。然而实测发现，对于三级BLAS，仅仅通过改变运算访存比性能提升空间还是有限的(仅为5％左右)，优化需要从运算本身考虑。

复数经典算法是针对两个复数之间的乘法，设复数A＝x+yi，B＝m+ni，C＝p+qi,计算新的C的值为A×B+C，计算x乘以m的值e(e←x×m)，计算y乘以n的值f(f←y×n)，计算p加e后减f的值，并将此值赋给p(p←p+e-f)，计算x加y的值k(k←x+y),计算m加n的值t(t←m+n)，计算q加k乘以t的积再减去e减去f的值，并将此值赋给q(q←q+k×t-e-f)，则得到新的C的值。这种算法可以将4乘4加个浮点操作转化成3乘7加。由于BLAS矩阵分块的特性，可以考虑将该算法用在块状矩阵的运算中，这样块内4×(nb)³+4×(nb)²个浮点运算就转化为3×(nb)³+7×(nb)²个，不考虑访cache的影响，当3×(nb)³+7×(nb)²<4×(nb)³+4×(nb)²即nb>3时，整体运算规模就会减小从而性能会有一定的提升，由于算法的改变，在分块大小不变的前提下，理论上性能提升比应该为

$\frac{n-b}{4\&times;\mathrm{nb}+4}\&times;100\%---\left(1\right)$

由公式(1)可见，矩阵分块大小nb值越大，性能提升越多，且比值接近25％。当nb＝48时提升比例接近23％时，其中每次计算会使用4×(nb)²个临时存储空间，该存储区可以在下一个分块时被循环反复使用，所以其空间开销可以忽略不计。

5、应用效果及应用推广

图2为采用本发明复数矩阵的优化方法在龙芯处理器上的优化效果仿真图。以常用数学库基础线性代数程序集(BLAS)中一个名为zGEMM的广义矩阵相乘的函数为例，图2中下面的曲线A是原来的zGEMM函数在龙芯上运行的运行速度，性能在700Mflops到800Mflops之间，图2中中间的曲线B是通过使用乘加指令、循环展开、矩阵分块、指令调度、寄存器轮转和数据预取等优化技术在龙芯处理器上得到的浮点吞吐率，zGEMM性能保持在1200Mflops到1300Mflops之间。图2中上面的曲线C是在中间的曲线B的基础上使用复数乘法优化和改变矩阵分块大小得到的运算速度，GEMM性能保持在1900Mflops到2100Mflops之间，它将中间的优化曲线B的浮点吞吐率提升了大约50％，是原zGEMM函数运算速度的2.6倍以上。

本发明复数矩阵的优化方法，是针对龙芯3A的四发射体系结构和特定的一级数据缓存大小进行的优化，由于计算了在龙芯体系结构下特定的展开粒度，对复数矩阵进行四乘四循环展开和针对龙芯的特定大小分块，将矩阵分块大小nb选取为一个小于52并使得24乘以矩阵分块大小nb的平方(nb)²小于龙芯处理器一级数据缓存64千字节(kb)的最大的值，得到在龙芯上最优的理想矩阵分块大小；并利用了数据存储的连续性和局部性对矩阵乘法中的矩阵进行合理划分与合并，来减少龙芯的一级数据缓存cache访问次数；对复数矩阵运算中的普通复数乘加运算利用了复数经典算法中两个复数之间的乘法算法，从而减小了运行规模，使得复数矩阵乘法运算在龙芯3A上的计算性能提高了大约50％，可使BLAS库在龙芯3A上的运算速度达到原来的2.5倍以上。

Claims (1)

1.一种在龙芯上的复数矩阵的优化方法，其特征在于操作步骤如下：

第一步：利用展开粒度公式K_AK_B/(K_A+K_B)＝f/g计算最优取数比，其中K_A和K_B分别为矩阵A和矩阵B的展开粒度，f/g为最优运算取数比，令两个矩阵展开粒度的乘积除以两个矩阵的展开粒度之和所得的值与最优取数比相同，在龙芯上调整两个矩阵的展开粒度K_A和K_B，使得在循环结构中保证矩阵A和矩阵B的运算取数比f/g为2：1，则所得的展开粒度即为在运算中一次循环需要展开的运算次数；

第二步：在矩阵运算时，对于矩阵内数据大于500的大规模矩阵计算，对每一分块进行运算时在循环中采用4×4的展开，每个展开的运算对应龙芯的一个发射槽，通过发射槽发射运算指令给龙芯运算器进行乘加运算；根据矩阵分块公式24×(nb)²≤NK，其中NK为龙芯的一级数据缓存cache的大小，将矩阵分块大小nb选取为一个小于52并使得24乘以矩阵分块大小nb的平方(nb)²小于龙芯处理器一级数据缓存64千字节的最大的值，得到在龙芯上最优的理想矩阵分块大小；

第三步：对于在复数矩阵运算中出现的普通复数乘加运算，利用复数经典算法中两个复数之间的乘法算法，设第一复数为A＝x+yi，第二复数为B＝m+ni，第三复数为C＝p+qi,计算第一复数A乘以第二复数B后再加第三复数C的值，将此值赋给第三复数作为第三复数的新值，算法流程如下：计算第一复数A的实部x乘以第二复数B的实部m的值为第一实数e←x×m，计算第一复数A的虚部y乘以第二复数B的虚部n的值为第二实数f←y×n，计算第三复数C的实部p加第一实数e后减去第二实数f的值，并将此值赋给第三复数C的实部p←p+e-f，计算第一复数A的实部x加第一复数A的虚部y的值为第三实数k←x+y,计算第二复数B的实部m加第二复数B的虚部n的值为第四实数t←m+n，计算第三复数C的虚部q加第三实数k乘以第四实数t的积再减去第一实数e和第二实数f的值，并将此值赋给第三实数k的虚部q←q+k×t-e-f，则得到第三复数C的新值，通过龙芯发射槽传给运算器参与分块矩阵的运算。