0 引言

    在图像视频处理和机器学习领域，矩阵运算规模达数千维甚至上兆维^[1]。矩阵运算，尤其是矩阵乘法O(n³)成为影响上述应用实时性的关键。研究低成本和低时间开销的大规模矩阵乘法求解方法具有极强的工程实用价值。

1 相关工作

    C=A×B(A、B和C矩阵均为M阶)的矩阵乘法的结果数据之间无关，故C阵的元素可同时计算。利用这一特性，多种改进矩阵乘法器方法被提出获得低时间复杂度，较为典型的有：STRASSEN V等人提出了利用张量代数实现的公式化矩阵乘法，降低矩阵乘的时间复杂度的Strassen算法^[2]；CANNON L E等人介绍了一种网格并行的Cannon算法^[3]，通过矩阵分块和地址循环的方法实现矩阵乘，最快执行时间为O(n)；此后，FOX G C等人提出一种基于超立方体结构的Fox算法^[4]；KUNG H T等人提出了一种阵列结构的脉动算法^[5]；沈俊忠、田翔等人提出了两种基于脉动的二维阵列实现方案^[6-7]。

    理论上采用上述矩阵阵列乘法器阵列规模越大，加速效果越显著，但是实际硬件实现中阵列规模受两个因素制约。

    首先，理论上脉动阵列的IO带宽应满足3fMB(这里f为工作频率，M为阵列规模，B为数据位宽，当f和B一定时，可记为O(3M)），受存储墙影响，M无法得到K级；

    其次，过高的M限制了阵列IO带宽利用率。设脉动阵列输入带宽利用率B_ave定义为：

    式(1)的函数曲线如图1所示，图1表明随阵列规模M的增长，B_ave值逐渐收敛为1/2。这就意味着单纯增加M值会造成近一半IO带宽的浪费，这又进一步加剧了IO带宽的压力。

    为克服经典脉动阵列结构的缺陷，KUMAR V K P等人提出一种IO带宽固定的链式乘法器^[8-9]，该乘法器中每个运算单元内存储容量固定，通过将矩阵A元素和矩阵B元素送到“低速通道”和“快速通道”来实现Block内和Block间的数据传递，这种数据组织形式将IO带宽限制为9fB，避免了随着矩阵规模增大带来的带宽问题，该结构另外的优点在于运算单元内部的存储容量固定，但容量和运算单元数量呈反比；对于大规模矩阵乘，容量越小意味着运算单元数量越多，造成运算单元数量的浪费，并且由于线性乘法器仅有单一的数据输入端，这会造成从数据输入到所有运算单元都处于工作状态的延时更大，降低了大规模矩阵运算效率。

    ALOGEELY M A等人提出的新型脉动阵列结构^[10]改进了KUMAR V K P等人提出的乘法器^[8-9]，通过对每个Block内的运算单元数量做了裁剪，构成梯度式增长的链式乘法器，这种结构解决了运算启动的延时问题，但数据组织较为复杂，不具备良好的拓展性。

    针对上述问题，本文提出了一种低IO带宽矩阵阵列乘法器——链式矩阵乘法器，主要特点是用列向量乘行向量和分时累加操作替代行向量乘列向量运算，使得乘法器内所有运算单元都参与每条向量计算过程，充分发挥了数据和运算单元可复用特性。本文给出了所设计链式矩阵乘法器的工作原理和对应硬件电路架构，并在FPGA芯片上完成乘法器硬件实现和性能测试。多种规模的矩阵乘法实测结果证实，本设计在极低IO带宽下达到了经典阵列乘法器计算性能。

2 链式乘法器工作原理

    矩阵乘法运算C=A×B，(A、B和C均为M维矩阵)的伪码如图2所示。

    令C[0][i][j]=0；1≤i≤M；1≤j≤M，设M=3，可得如图3左侧所示的三阶矩阵乘的DAG，图中每个顶点代表一个乘累加运算，无括号坐标表征图2中数据移动路线，该DAG的硬件实现如图3右侧所示。

    当IO带宽满足，规模为n×n的脉动乘法器可得计算时间下界O(n)。但受IO带宽制约，典型脉动阵列规模有限，需要将大维度矩阵分解若干子阵分别处理，执行复杂的数据缓冲或访存操作才能得到理想的加速效果。

    将图3左侧图沿d₁=(0，0，1)方向投影可得图4。

    图2中沿着i坐标轴正方向存在三个与坐标轴jk平面平行的平面，第二个平面的输入是行向量a₂和矩阵B，第三个平面输入是行向量a₃和矩阵B。故其运算过程可描述为如下公式：

    式(2)中，矩阵A、B分别以列主式和行主方式展开乘法运算。由此可构造一种无源数据缓存的3个数据通道(分别对应源矩阵A和B，结果矩阵C)的链式乘法器，其计算过程的伪代码如图5所示，其运算路线如图2中带括号的坐标。

    综上，本文提出了一种IO带宽恒为3fB的链式乘法器，矩阵A、B中元素分别按列和行输入乘法器内运算单元完成确定操作。该链式乘法器支持在线配置，能够根据待处理矩阵规模被配置多链结构，适应不同的带宽条件。

3 链式矩阵乘器硬件实现

3.1 链式矩阵乘法器架构设计

    由图4得到行向量a₁与矩阵B相乘的数据组织形式如图6(a)所示。图中大方框中是用来执行乘累加操作并缓存结果矩阵元素的2组中间值，进而输出最终结果。由此类推，行向量a₂和a₃与矩阵B相乘的数据组织形式如图6(b)和图6(c)所示。

    由图6可知，M阶(此时M=3)矩阵A在与M阶矩阵B进行矩阵乘时，矩阵A中的行向量的每个元素都要保持M个周期才能向对应乘累加器中输入下一个元素，即输入带宽利用率仅为矩阵A输入总带宽的1/M；矩阵B按行向每个乘累加器输入相同的元素，其带宽利用率也仅为矩阵B输入总带宽的1/M。观察矩阵A的行向量的输入规律，在行向量的每个元素的保持周期内，将矩阵A中的元素按列进行输入(例如，在矩阵A第一个行向量的第一个元素a₁₁的保持周期内，将矩阵A中的元素a₁₁、a₂₁、a₃₁依次输入到每一个乘累加器中)，可以得到三阶矩阵列乘行的数据组织形式如图7所示。图7中矩阵A按列依次进行输入，矩阵B按行依次进行输入，但要使得矩阵B输入端的带宽得到充分的利用，矩阵B需要在每个相邻乘累加器相差1cycle，故而可以将矩阵B也按照行输入到每个乘累加器，通过使用寄存器将矩阵B的输入延迟一个周期输入到下一个乘累加器中，可以得到图8的数据组织形式。

3.2 链式乘法器PE设计

    如图8所示，为了将链式乘法器的输出带宽利用率达到最高，本文设计了如图9所示PE的结构，包含一个浮点数乘法器和一个浮点数加法器，一对用于脉动输入的寄存器，以及一组深度为m的乒乓存储器，分别用于缓存加法器结果和输出矩阵运算结果。

    当源数据srcA和srcB进入PE时，首先被送入浮点数乘法器中进行乘法运算，乘法结果会被送入浮点数加法器，与运算存储器(用于缓存加法器结果的存储器)中顶层(dstA0)数据相加，相加的结果进入运算存储器底层(dstA(m-1))；加法器每进行一次运算，所有数据被刷新一次，即向前移动一个地址，直到处于顶层(dstA0)；反复执行上述操作，直到运算完毕。

    当前矩阵运算结束后，需要将结果输出，此时存储器执行乒乓操作，将运算存储器和原先用于输出矩阵运算结果的传输存储器互相切换。

    这种结构将累加操作解耦，只将乘法和加法操作结合，使得PE更加独立，同时避免了硬件实现中累加器的流水线级数造成的结果数据延迟效应。

3.3 链式乘法器硬件设计

    以M阶矩阵为例，链式乘法器需要M个PE。在进行C=A×B的矩阵乘法时，输入矩阵A和B分别按列和按行输入到链式乘法器结构中。图10所示为链式乘法器结构中M=3的具体结构，矩阵B的数据元素每周期依次进入PE₀的1端口，而矩阵A中的元素每周期依次对各个PE的0端口轮转输入。

    链式乘法器的工作流程如图11所示，具体对应为一个3阶方阵相乘的流程。可以看到，在初始M(M=3)个周期的启动时间过后，所有PE处于满载状态。对于PE₀，直到第8个周期才输出结果c11，c11是由第1、4、7周期的元素乘加得到的，三个阶段的结果分别对应图中c11₁、c11₂、c11₃，最终求和得到的结果为矩阵A的第1行向量与矩阵B的第1列向量的乘累加结果。图中灰色部分表示新的矩阵A′和B′输入并相乘，此时的结果传输时间被隐藏。

    根据这种矩阵乘法的结构特性，可以实现多个矩阵的连乘操作。由于没有输入缓冲，一组矩阵完成计算并切换到下一矩阵时没有存储上的时间消耗，结果存储器采用乒乓操作，一部分用于参与当前矩阵的实时运算，另外一部分存储了上次矩阵乘结果并回传，隐藏了大部分数据回写时间。

    在面对大规模矩阵乘的时候，由于硬件资源受限，难以一次性完成整个矩阵的运算，这时候必须将矩阵分成多个小块，对各个子块以及块间进行乘加运算，链式乘法器能够很好地处理分块矩阵乘。

4 链式乘法器性能分析

    下面以在FPGA上实现的链式乘法器来对上述设计的性能进行分析，在本文中，设计在XC7V2000T芯片上进行验证，实验中采取单精度浮点数作为数据元素，经过实验，片内最大可集成800个PE。本设计中所使用的FPGA开发环境和仿真环境为Xilinx Vivado 2016.3及Mentor Graphics Modelsim SE-64 10.2c。

4.1 峰值计算性能

    理想情况下，链式乘法器工作时所有PE均满载，每个PE在单时钟周期内消耗两个浮点数，得到一个结果数据。实际上，由于矩阵类型、流水线延迟、存储策略加上启动延迟(数据从输入开始到所有PE工作)等因素的影响，在某些周期内仍会存在部分PE处于空闲的情况，设PERF_peak表示运算的峰值性能，N表示PE的个数，f表示工作频率，S表示整个大规模矩阵总运算次数，则有式(3)表示如下：

    本文主要围绕大规模矩阵储乘法进行硬件设计，通过探究N、S、f等因素对运算器的性能影响，进行了以下几组实验，实验以200 MHz时钟作为工作时钟，以2片DDR3作为主存，存储源数据和结果数据，2片DDR3峰值带宽为204.8 Gb/s，考虑到多路并行输入输出时的带宽限制问题，在实验中矩阵元素采用32位浮点数表示，所以单条链式乘法器的峰值带宽为18.75 Gb/s(2组源数据，一组结果数据)。

    在xc7v2000tflg1925-1 FPGA上实现该设计，FPGA内部资源使用情况如表1所示，根据布线后仿真的结果，该矩阵乘法器在未做优化的情况下工作频率可达到100 MHz。由此得出，该设计的峰值单精度浮点计算性能可以达到156.25 GFLOPS。

    当PE资源不足以支持一次性运算整个矩阵，需要分块运算，表2中阴影部分表示运算不需分块，非阴影部分表示需要分块运算。可以看出，随矩阵规模增大，不同链数对运算性能的影响逐渐减小，图中，在处理128×128及更大规模的矩阵乘时，所有链式乘法器运算周期基本持平，这是因为在处理这些矩阵分块运算时，所有PE均处于满负荷状态，而矩阵分块仅仅改变的是同一数据的读取次数，对总的运算量并不影响。

4.2 不同链数对运算的影响

    分析表2中所列的数据，可以得到，当运算规模小的时候，链式乘法器的整体运算时间相较于数据运算量来说较大，这是因为数据回传时间对于总运算周期占比较大，随着矩阵规模的增大，数据回传时间占总周期比重越来越小，计算效率逐渐提高。

    由于表2中数据范围相差较大，为了直观地得到各组乘法器的性能对比，对每组链的运算周期与理想脉动结构的周期做比值,得到结果如图12所示。

    表2中同样列举了规模为8×8的脉动阵列的运算周期，200 MHz工作频率下8×8脉动阵列的峰值带宽为150 Gb/s，而本文设计的单链乘法器峰值带宽仅18.75 Gb/s，由图12可知，单链乘法器在处理小规模矩阵时与脉动相比性能有所欠缺，但在处理大规模矩阵乘时两者表现相当，而单链乘法器在带宽方面优势非常明显。

    设PPB表示链式结构单位带宽的运算能力，并以脉动阵列作为衡量标准，PPB为脉动阵列所有带宽的运算能力与链式结构的运算能力的归一化比值，满足如下关系：

    式(4)中，链数L₁表示脉动阵列规模(本实验中L₁=8)，L₂表示链式乘法器的链数，T_systolic和T_chain分别表示脉动和链式结构的运算周期，将表2中数据带入式(2)，得到链式乘法器的PPB如图13所示。

    当运算小规模矩阵时，链式乘法器单位带宽运算能力较低，然而当计算大规模矩阵乘时，由于链式乘法器能够很好地处理矩阵分块造成的运算效率问题，因此单位带宽的运算能力逐渐提高，并逐渐趋近于L₁与L₂的比值。例如本实验中，单链模式下，PPB最终趋近于8，表示链式乘法器单位带宽的运算能力与脉动阵列的8带宽的运算能力相当。由图13可知，本文所设计的链式乘法器在单链(L=1)下工作性能最佳，同比与其他模式下的运算器，单链乘法器在带宽方面有着明显的优势。

4.3 不同数量运算单元对运算的影响

    不同规模的链式乘法器对于矩阵乘法也有不同的加速效果，本文分别以PE集成度为16、32、64、128的单链乘法器为实验对象，计算不同规模的矩阵乘，见表3。

    表3中的数字表示总运算周期，阴影部分表示不分块运算结果，非阴影部分表示分块运算结果。可以看出，在处理小规模矩阵时，所有运算器性能相当，这是因为PE资源足够支撑一次性运算；在处理大规模矩阵时，PE数量成为限制运算器性能的主要因素，在处理128×128及更大规模的矩阵乘时，运算周期随PE数量呈递减趋势。

    分块运算时，受PE数量N限制，每次运算量为N×N，若将矩阵分割为K×K个N×N的矩阵块，则总的运算周期近似N²×K³。设其他结构的乘法器所分块得到的子阵规模为N′×N′，其分块得到的子阵数量为K′×K′，可以得到，相同PE数量下，单链结构所对应的N值大于N′，K小于K′，所以总运算周期最优。

    综上所述，本文提出的链式乘法器可以处理不同规模的矩阵乘法，并且运算性能卓越。

5 结论

    本文设计了一种链式并行浮点数矩阵乘法器，并在Xilinx的XC7V2000T芯片进行了原型验证，通过实验分析，该矩阵乘法器优点在于对IO带宽的要求非常低，结构灵活，能够适用于不同类型的矩阵乘法。该设计主要面向大规模矩阵的分块运算，由于数据采取非缓冲的组织形式，运算器能够实现数据流入的同时开始计算，数据流入完毕结束运算，极大地提高了整体运算的吞吐率；由于PE相对独立，支持根据不同的矩阵规模根据配置信息在线调整PE的使用情况，满足了灵活性与通用性的要求，同时具备了良好的拓展性，链式乘法器在运算大规模矩阵乘法时表现突出。

参考文献

[1] 冯子勇.基于深度学习的图像特征学习和分类方法的研究及应用[D].广州：华南理工大学，2016.

[2] STRASSEN V.Gaussian elimination is not optimal[J].Numerische Mathematik，1969，13(4)：354-356.

[3] CANNON L E.A cellular computer implement the Kalman filter algorithm[D].Bozeman US：Montana State University，1969.

[4] FOX G C，OTTO S W，HEY A J G.Matrix algorithm on a hypercube I: matrix multiplication[J].Parallel Computing，1987，4(1)：17-31.

[5] KUNG H T，LEISERSON C E.Systolic arrays(for VLSI)[J].Proceeding Sparse Matrix Conference，1978：256-282.

[6] 沈俊忠，肖涛，乔寓然，等.一种支持优化分块策略的矩阵乘加速器设计[J].计算机工程与科学，2016，38(9)：1748-1754.

[7] 田翔，周凡，陈耀武，等.基于FPGA的实时双精度浮点矩阵乘法器设计[J].浙江大学学报(工学版)，2008，42(9)：1611-1615.

[8] KUMAR V K P，TSAI Y C.On synthesizing optimal family of linear systolic arrays for matrix multiplication[J].IEEE Transactions on Computers，1991，40(6)：770-774.

[9] KUMAR V K P，TSAI Y C.On Mapping algorithms to linear and fault-tolerant systolic arrays[J].IEEE Transactions on Computers，1989，38(3)：470-478.

[10] ALOGEELY M A，Al-Turaigi M A，ALSHEBEILI S A.A new approach for the design of linear systolic arrays for computing third-order cumulants[J].Integration the Vlsi Journal，1997，24(1)：1-17.

作者信息:

宋宇鲲，郑强强，王泽中，张多利

(合肥工业大学 电子科学与应用物理学院，安徽 合肥230009)