一种支持可变分块的矩阵乘加速方法
一种支持可变分块的矩阵乘加速方法
【技术领域】
[0001] 本发明涉及嵌入式平台下的矩阵乘加速技术，具体涉及一种支持可变分块的矩阵 乘加速方法。
【背景技术】
[0002] 随着半导体制造工艺的发展以及集成电路技术的进步，单芯片上能够集成越来越 多的晶体管，使用可编程器件，特别是FPGA(FieldProgrammableGateArray)芯片进行设 计成为了目前构建嵌入式系统以及硬件加速平台的一种重要方式。当前FPGA芯片提供了 专用的算术模块、大量的逻辑资源和存储资源，以及外部存储器接口、网络接口和其它外围 接口，为构建高性能计算系统提供了条件，也使FPGA可重构计算系统成为加速科学计算的 一种重要选择。当前的FPGA芯片相对于DSP(DigitalSignalProcessor,数字信号处理 器）具有可编程的优势，同时能够并行处理海量数据，既具有通用处理器的灵活性，又具有 ASIC(ApplicationSpecificIntegratedCircuit,专用集成电路）的高性能，在嵌入式计 算领域备受青睐。
[0003] 浮点矩阵乘法是数字信号的基本算法，同时也是许多科学计算方法的基本运算。 在数字图像处理，计算机视觉的快速处理以及工业实时控制等领域都被广泛应用。但由于 实际应用中，浮点矩阵规模通常较大，矩阵乘算法本身复杂度较高、处理效率较低，成为限 制系统性能提升的瓶颈所在，因此为此类应用设计高性能的硬件结构是当前FPGA结构设 计的研宄热点。
[0004] 尽管近年来提出了很多基于FPGA的矩阵乘加速器设计，但是都缺乏对非均匀的 大规模矩阵加速的讨论和支持，这种大规模矩阵的特征是行列数相差很大（> =10倍），并 且广泛存在于很多现代应用领域中，如图像处理，深度学习等。在这些应用领域之中，矩阵 乘占据了计算量的主要部分。由于单片FPGA芯片上存储资源和计算资源十分有限，加速大 规模矩阵乘时往往需要对矩阵进行分块。对于链式结构的矩阵乘加速器，虽然对大多数大 规模矩阵的加速效果十分明显，但对于加速非均匀矩阵时的计算效率却很低，其主要原因 就是这类加速器往往只支持固定的分块，也就是说，分块大小与矩阵链长（矩阵链中处理 单元个数）相等或者是其倍数。当分块大小与矩阵加速器链长不匹配时，加速器计算效率 会发生明显的下降。据我们所知，至今还没有公开文献涉及支持可变分块的矩阵乘法器设 计，也没有经典分块算法基础上的分块优化问题方面的相关研宄，因此如何选择最优的分 块，使得矩阵乘加速器获得更高的计算效率，以更好地适应现代应用的加速需求，已经成为 一项亟待解决的关键技术问题。
【发明内容】
[0005] 本发明要解决的技术问题是：针对现有技术的上述技术问题，提供一种支持可变 分块，可根据分块大小调整使用的处理单元数目，加速非均匀矩阵乘运算加速效率高的支 持可变分块的矩阵乘加速方法。
[0006] 为了解决上述技术问题，本发明采用的技术方案为：
[0007] -种支持可变分块的矩阵乘加速方法，步骤包括：
[0008] 1)输入矩阵乘运算所需的M*N的矩阵A和N*R的矩阵B;
[0009] 2)根据矩阵A和矩阵B的规模确定子块大小Si，将矩阵A以规模为Si*N的子块为 单位进行按行分块，将矩阵B以规模为N*Si的子块为单位进行按列分块，使得矩阵乘运算 等同为多个子块乘运算；
[0010] 3)为每一个子块乘运算所需数据生成一个DMA描述符，将所有子块乘运算的DMA描述符构建DMA描述符链表并存入主存；
[0011] 4)针对每一个子块乘运算，通过矩阵乘加速器的DMA从主存的DMA描述符链表读 取子块乘运算所需数据，然后通过矩阵乘加速器中至少一条处理单元链中的前Sif处理单 元进行子块乘运算，并将各个子块乘运算的结果通过矩阵乘加速器的DMA分别写回主存。
[0012] 优选地，所述步骤2)中根据矩阵A和矩阵B的规模确定子块大小Si的详细步骤 包括：
[0013] 2. 1)输入矩阵A和矩阵B的矩阵规模，所述矩阵规模包括矩阵A的行数M、矩阵A 的列数N、矩阵B的列数R三者的值；
[0014] 2. 2)根据矩阵规模获取满足式⑴所示约束的分块大小集合；
[0016] 式（1)中，3,表示矩阵A被按行划分的行数，S」表示矩阵B被按列划分的列数， M表示矩阵A的行数，N表示矩阵A的列数，BW表示矩阵乘加速器中单条处理单元链的 访存带宽，max以，Sj表示从SJP\中取较大值，F表示访存频率；P表示矩阵乘加速 器中单条处理单元链的处理单元数量，Stageadd表示矩阵乘加速器中的加法器流水级数， maxm-SiXi^R-SjXn}表示从M-SiXm和R-SjXn中取较大值，R表示矩阵B的列数，m表示 矩阵A的行数M除以矩阵A被按行划分的行数结果的向下取整值，n表示矩阵B的列 数R除以矩阵B被按列划分的列数&的结果的向下取整值；
[0017] 2. 3)将矩阵A视为包含整数个规模为S#N的子块的子矩阵①和剩余的不规则的 子矩阵②、将矩阵B视为包含整数个规模为N*Si的子块的子矩阵③和剩余的不规则的子矩 阵④，将矩阵A和矩阵B的矩阵乘运算视为子矩阵①~子矩阵④四者中的两两相乘，建立式 (2)所示的评估函数；
[0018]f(S^ Sj)=T1； 3(S^ Sj)+T1； 4(S^ Sj)+T2； 3(S^ Sj)+T2； 4(S^ Sj)(2)
[0019] 式⑵中，f(Si，Sp表示评估函数的值，5,表示矩阵A被按行划分的行数，5」表 示矩阵B被按列划分的列数，USpSp表示子矩阵①乘子矩阵③的计算时钟节拍数， USi，Sj表示子矩阵①乘子矩阵④的计算时钟节拍数，T2,3(Si，Sj表示子矩阵②乘子矩阵 ③的计算时钟节拍数，T2,4(Si，SP表示子矩阵②乘子矩阵④的计算时钟节拍数；
[0020] 2. 4)使用MATLAB数学工具对式（2)所示的评估函数生成满足式（1)约束的可视 化图像，以矩阵A被按行划分的行数Si和矩阵B被按列划分的列数Sd目等为前提，通过所 述可视化图像确定评估函数定义域内的全局最小值，得到分块大小集合内最优的分块大小 si〇
[0021] 优选地，所述步骤3)的详细步骤包括：
[0022] 3. 1)将矩阵A以规模为SjN的子块为单位按行分块，得到包含整数个规模为SjN 的子块的子矩阵①，如果仍有剩余的不规则的子块，则将剩余的不规则的子块作为子矩阵 ②；将矩阵B以规模为N*Sd^子块为单位按行分块，得到包含整数个规模为N*S,的子块的 子矩阵③，如果仍有剩余的不规则的子块，则将剩余的不规则的子块作为子矩阵④；将矩阵 A和矩阵B的矩阵乘运算视为得到的所有子矩阵的两两相乘，根据得到的子矩阵数量建立 使用指针相连接的DMA描述符链表，使得每一对相乘的子矩阵对应一个DMA描述符链表；
[0023] 3. 2)选择一个子块乘运算作为当前子块乘运算；
[0024] 3. 3)获取当前子块乘运算在矩阵A中对应子块数据的首地址、传输长度及传输步 长，获取当前子块乘运算在矩阵B中对应子块数据的首地址、传输长度及传输步长，将当前 子块乘运算在矩阵A中对应子块数据的首地址、传输长度及传输步长和当前子块乘运算在 矩阵B中对应子块数据的首地址、传输长度及传输步长封装生成一个DMA描述符，根据当前 子块乘运算所属的一对相乘的子矩阵确定对应的DMA描述符链表，并将该DMA描述符写入 对应的DMA描述符链表中；
[0025] 3. 4)判断是否所有子块乘运算已经完成处理，如果尚未完成所有子块乘运算的处 理，则选择下一个子块乘运算作为当前子块乘运算，跳转执行步骤3. 3);否则，如果已经完 成所有子块乘运算的处理则跳转执行步骤4)。
[0026] 优选地，所述步骤4)的详细步骤包括：
[0027] 4. 1)CPU将DMA描述符链表中第一个DMA描述符的首地址配置给矩阵乘加速器的 DMA，矩阵乘加速器的DMA开始根据首地址读取第一个DMA描述符作为当前描述符；
[0028] 4. 2)矩阵乘加速器的DMA解析当前描述符，得到当前描述符对应子块乘运算所需 的一对子块在外存中的地址及传输长度，并根据预设的传输步长和当前描述符中携带的一 对子块的地址及传输长度，以分时的方式交换读取当前描述符对应子块乘运算所需数据， 并将输入数据存入FIFO缓存；
[0029] 4. 3)基于FIFO缓存中当前描述符对应子块乘运算所需数据，通过矩阵乘加速器 中至少一条处理单元链中的前Si个处理单元进行当前描述符对应子块乘运算；
[0030] 4. 4)将当前描述符对应子块乘运算的结果通过矩阵乘加速器的DMA分别写回主 存；
[0031] 4. 5)判断DMA描述符链表中的所有DMA描述符是否已经处理完毕，如果尚未处理 完毕，则从DMA描述符链表中选择下一个DMA描述符作为当前描述符，跳转执行步骤4. 2); 如果已经处理完毕，则结束并退出。
[0032] 优选地，所述步骤4. 3)的详细步骤包括：
[0033] 4. 3. 1)在矩阵乘加速器中确定一条用于当前描述符对应子块乘运算的处理单元 链，所述处理单元链至少包含Si个处理单元；所述矩阵乘加速器中设有用于控制各个处理 单元工作状态的状态机，所述状态机包括预取阶段、预取及计算阶段、提交阶段共三个状 态；初始化变量k为0,控制状态机进入预取状态；
[0034] 4. 3. 2)在预取阶段，针对当前描述符对应子块乘运算，矩阵乘加速器的DMA将属 于矩阵A的子块中的第k列数据附带上Sif处理单元的编号信息 后发出给处理单元链，处 理单元链的前Si个处理单元根据编号信息预取编号对应的数据，控制状态机进入预取及计 算阶段；
[0035] 4. 3. 3)在预取及计算阶段，矩阵乘加速器的DMA将属于矩阵B的子块中第k行数 据以数据流的形式依次通过处理单元链的前Si个处理单元，且前Sif处理单元每一拍将接 收到的数据和属于矩阵A的子块中的第k列数据部分进行乘法运算，并将乘法运算结果和 片上存储中存储的上一次累加结果进行累加，将累加结果写回片上存储作为下一次累加的 输入；同时，处理单元链的前Sif处理单元分别根据编号信息预取属于矩阵A的子块中的 第k+1列数据；
[0036] 4. 3. 4)判断变量k的值是否等于N减1，如果变量k的值等于N减1，则控制状态 机进入提交阶段，跳转执行步骤4. 3. 5);否则，将初始化变量k加1，控制状态机进入预取及 计算阶段，跳转执行步骤4. 3. 3);
[0037] 4. 3. 5)在提交阶段下，处理单元链的前Si个处理单元分别将最后一次累加的结果 递给矩阵乘加速器的DMA，最终由矩阵乘加速器的DMA将结果写入主存中的指定区域。
[0038] 本发明支持可变分块的矩阵乘加速方法具有下述优点：本发明通过确定子块大小 Sp将矩阵A以规模为S#N的子块为单位进行按行分块，将矩阵B以规模为N*Si的子块为单 位进行按列分块，使得矩阵乘运算等同为多个子块乘运算，为每一个子块乘运算所需数据 生成一个DMA描述符，将所有子块乘运算的DMA描述符构建DMA描述符链表并存入主存，针 对每一个子块乘运算，通过矩阵乘加速器的DMA从主存的DMA描述符链表读取子块乘运算 所需数据，然后通过矩阵乘加速器中至少一条处理单元链中的前Sif处理单元进行子块乘 运算，并将各个子块乘运算的结果通过矩阵乘加速器的DMA分别写回主存，增加了对可变 分块的支持，能够根据具体的分块大小调整使用的处理单元数目，通过确定子块大小3 1来 基于经典分块算法上的优化分块策略，当加速非均匀矩阵时，通过计算并应用最优分块能 够使得矩阵加速器具有很高的计算效率，具有支持可变分块，可根据分块大小调整使用的 处理单元数目，加速非均匀矩阵乘运算加速效率高的优点。
【附图说明】
[0039] 图1为本发明实施例的基本流程示意图。
[0040] 图2为本发明实施例将进行矩阵乘运算等同为多个子块乘运算的原理示意图。 [0041] 图3为本发明实施例将进行矩阵乘运算划分为子矩阵相乘的原理示意图。
[0042] 图4为本发明实施例中使用MATLAB数学工具生成的评估函数可视化图像。
[0043] 图5为本发明实施例中生成的一个DMA描述符链表的结构不意图。
[0044] 图6为本发明实施例中矩阵加速器的DMA的工作流程示意图。
[0045] 图7为本发明实施例中加速器系统的框架结构示意图。
[0046] 图8为本发明实施例中处理单元的数据流向结构示意图。
[0047] 图9为本发明实施例和采用经典分块算法的性能对比图。
【具体实施方式】
[0048] 如图1所示，本实施例支持可变分块的矩阵乘加速方法的步骤包括：
[0049] 1)输入矩阵乘运算所需的M*N的矩阵A和N*R的矩阵B;
[0050] 2)根据矩阵A和矩阵B的规模确定子块大小Si，将矩阵A以规模为Si*N的子块为 单位进行按行分块，将矩阵B以规模为N*Si的子块为单位进行按列分块，使得矩阵乘运算 等同为多个子块乘运算；
[0051] 3)为每一个子块乘运算所需数据生成一个DMA描述符，将所有子块乘运算的DMA 描述符构建DMA描述符链表并存入主存；
[0052] 4)针对每一个子块乘运算，通过矩阵乘加速器的DMA从主存的DMA描述符链表读 取子块乘运算所需数据，然后通过矩阵乘加速器中至少一条处理单元链中的前Sif处理单 元进行子块乘运算，并将各个子块乘运算的结果通过矩阵乘加速器的DMA分别写回主存。
[0053] 如图2所示，对于M*N的矩阵A和N*R的矩阵B的矩阵乘运算而言，其矩阵乘运算 得到的矩阵C的规模为M*R。因此，将矩阵A和N*R的矩阵B的矩阵乘运算划分为多个小 矩阵相乘，其中矩阵A按行划分为多个S#N子块，矩阵B按列划分为多个N*Sj子块，每对小 子块相乘就得到一个SjSj勺矩阵，也就是最后的结果C矩阵的一个子块。参见图2,现考 虑子块相乘的细节，可把A的S#N子块看成是由N个维度为列向量组成的矩阵，B的 N*Sj子块是由N个维度为S」的行向量组成的矩阵，记第k个列向量和行向量分别为Vk和 Uk(k〈 =N)，把Uk和Vk相乘的结果进行累加（N次）就得到矩阵C的SjSj子块。假设矩阵 乘加速器有Si个处理单元，每个处理单元每次缓冲Vk中的一个操作数，Uk中的所有操作数 依次进入每个处理单元，处理单元以流水线的形式高效地进行乘加运算，并且将加法结果 存入临时存储区（下一次加法运算要从存储区读取这一次的结果进行累加），这样一来，每 个处理单元并行地计算C的某一行的结果，这就是本实施例算法的核心思想。事实上当Si =SJt，矩阵乘加速器进行阶段同步所导致的空拍时最少，这种情况下矩阵乘加速器的性 能应该是最优的，因此本实施例仅考虑Si=Sj勺情况。
[0054] 本实施例中，步骤2)中根据矩阵A和矩阵B的规模确定子块大小Si的详细步骤 包括：
[0055] 2. 1)输入矩阵A和矩阵B的矩阵规模，矩阵规模包括矩阵A的行数M、矩阵A的列 数N、矩阵B的列数R三者的值；
[0056] 2. 2)根据矩阵规模获取满足式⑴所示约束的分块大小集合；
[0058] 式（1)中，31表示矩阵A被按行划分的行数，5』表示矩阵B被按列划分的列数， M表示矩阵A的行数，N表示矩阵A的列数，BW表示矩阵乘加速器中单条处理单元链的 访存带宽，max以，Sj表示从SJP\中取较大值，F表示访存频率；P表示矩阵乘加速 器中单条处理单元链的处理单元数量，Stageadd表示矩阵乘加速器中的加法器流水级数， maxm-SiXi^R-SjXn}表示从M-SiXm和R-SjXn中取较大值，R表示矩阵B的列数，m表示 矩阵A的行数M除以矩阵A被按行划分的行数结果的向下取整值，n表示矩阵B的列 数R除以矩阵B被按列划分的列数&的结果的向下取整值；
[0059] 2. 3)将矩阵A视为包含整数个规模为S#N的子块的子矩阵①和剩余的不规则的 子矩阵②、将矩阵B视为包含整数个规模为N*Si的子块的子矩阵③和剩余的不规则的子矩 阵④，将矩阵A和矩阵B的矩阵乘运算视为子矩阵①~子矩阵④四者中的两两相乘，建立式 (2)所示的评估函数；
[0060] f(S^Sj) =T1； 3 (S^Sj) +T1； 4 (S^Sj) +T2； 3 (S^Sj) +T2； 4 (S^Sj) (2)
[0061] 式（2)中，表示评估函数的值，51表示矩阵A被按行划分的行数，Sj-表 示矩阵B被按列划分的列数，USpSp表示子矩阵①乘子矩阵③的计算时钟节拍数， USi，Sj表示子矩阵①乘子矩阵④的计算时钟节拍数，T2,3(Si，Sj表示子矩阵②乘子矩阵 ③ 的计算时钟节拍数，T2,4(Si，SP表示子矩阵②乘子矩阵④的计算时钟节拍数；
[0062] 2. 4)使用MATLAB数学工具对式（2)所示的评估函数生成满足式（1)约束的可视 化图像，以矩阵A被按行划分的行数Si和矩阵B被按列划分的列数Sd目等为前提，通过可 视化图像确定评估函数定义域内的全局最小值，得到分块大小集合内最优的分块大小Sp
[0063] 如图3所示，将矩阵A和矩阵B的矩阵乘运算视为子矩阵①~子矩阵④四者中的 两两相乘后，子矩阵①包含整数个规模为SjN的子块，子矩阵③包含整数个规模为N*SJ9 子块，剩余的不规则的子矩阵②包含的为不规则额的子块（行数小于SJ，剩余的不规则的 子矩阵④包含的为不规则的子块（列数小于SJ，子矩阵①~子矩阵④四者中的两两相乘 即为子矩阵①乘子矩阵③、子矩阵①乘子矩阵④、子矩阵②乘子矩阵③、子矩阵②乘子矩阵 ④ 。本实施例中，由于子块乘运算要求在开始的时候预取Si个矩阵A的数据（子块第一 列），这个过程需要Si个节拍，每个PE计算SJSj需要花费maxtSpSjlXN个节拍。由于 矩阵乘加速器中的处理单元组织成链式结构，因此最后一个计算的处理单元相对于第一个 处理单元有Si个节拍的延迟，即整个处理单元流水线的建立时间。计算S，5」子块的结 果需要从主存读取（SJSPXN个数据，该过程需要时间（Si+SpXN/BWJW是单条链的访存 带宽。根据所有处理单元总计算时间大于取数时间，有以下不等式（1-1)成立。
[0064] (Sj+Sj)XN/BW^ (2XSj+max{S^Sj}XN)/F(1-1)
[0065] 式（1-1)中，3,表示矩阵A被按行划分的行数，S」表示矩阵B被按列划分的列数， N表示矩阵A的列数，BW表示矩阵乘加速器中单条处理单元链的访存带宽，max以，Sj表示 从SJP\中取较大值，F表示访存频率。
[0066] 此外，矩阵A被按行划分的行数Sp矩阵B被按列划分的列数&还有以下约束：首 先，根据算法，矩阵A被按行划分的行数Si的值不能大于矩阵乘加速器中一条处理单元链 中处理单元的数目；其次，考虑到流水线中的数据冲突，也就是当加法结果尚未写入片上存 储器时，加法器又需要从存储器读取该结果时，加法器就会读到旧值。因此，必须同时满足 式（1-2)所示的几个约束条件。
[0068] 式（1-2)中，51表示矩阵A被按行划分的行数，P表示矩阵乘加速器中单条处理单 元链的 处理单元数量，max以，S」}表示从SjPS」中取较大值，Stageadd表示矩阵乘加速器 中的加法器流水级数，11^復-5'111，1^，11}表示从1-5'111和1^，11中取较大值。本实 施例中，Stageadd表示矩阵乘加速器中的加法器流水级数具体值为11。结合前述的式（1-1) 和式（1-2)，即可推导得出本实施例前述的式（1)所示约束的分块大小集合。
[0069] 假设?q」，" =[RJ(即m表示矩阵A的行数M除以矩阵A被按行划分的行 数结果的向下取整值，n表示矩阵B的列数R除以矩阵B被按列划分的列数\的结果 的向下取整值），则有子矩阵①乘子矩阵③的总计算时钟节拍数USpSj如式（2-1)所 示，子矩阵①乘子矩阵④的总计算时钟节拍数USpSj如式（2-2)所示，子矩阵②乘子矩 阵③的总计算时钟节拍数T2,3(Si，Sp如式（2-3)所示，子矩阵②乘子矩阵④的总计算时钟 节拍数如式（2-4)所示，
[0070] 1\，3以，Sj) =mXnX(Si+maxtSi，SjXN) (2-1)
[0071] Sj) =k2XmX(maxtSj,R-SjXn}XN) (2-2)
[0072] T2； 3 (S^Sj) =k:XnX(max{M-SjXm,Sj}XN) (2-3)
[0073] T2； 4 (S^Sj) =k:Xk2X(max{M-SjXm,R-SjXn}XN) (2-4)
[0074] 式（2-1)~（2-4)中，maxtSi，S」}表示从SJPS」中取较大值，max{Si，R-SjXn} 表示从SJPR-S」Xn中取较大值，max{M-SiXm，S」}表示从M-SiXm和Sj中取较大值， maxm-SiXi^R-SjXn}表示从M-SiXm和R-SjXn中取较大值，m表示矩阵A的行数M除以 矩阵A被按行划分的行数结果的向下取整值，n表示矩阵B的列数R除以矩阵B被按 列划分的列数h的结果的向下取整值，S,表示矩阵A被按行划分的行数，S^表示矩阵B被 按列划分的列数，M表示矩阵A的行数，N表示矩阵A的列数，R表示矩阵B的列数。
[0075] 式（2-2)~（2-4)中，变量&和1^2的值由式（1-2)和（1-3)决定。
[0078] 式（1-2)和（1-3)中各个参数的定义可参见式（2-1)~（2-4)，在此不再赘述。
[0079] 参见式（2-1)~（2-4)可知，评估函数（2)的形式十分复杂，本实施例中仅考虑Si =S#情况。因此，根据变量kJPk2的值可能有下述情形：
[0080] 情形 1 %= 0,k2= 0):
[0081] 在这种情况下，m=M/Si，n=R/Sj，式（2)的评估函数可简化为式（2-5)。
[0083] 式（2-5)中，\3以，SJ表示子矩阵①乘子矩阵③的总计算时钟节拍数，M表示矩 阵A的行数，N表示矩阵A的列数，R表示矩阵B的列数，5,表示矩阵A被按行划分的行数。
[0084] 情形 2〇^=0,k2=l):
[0085] 在这种情况下，m=M/Sp 式⑵的评估函数可简化为式（2- 6)。
[0086]
[0087] 式（2-6)中，表示子矩阵①乘子矩阵③的总计算时钟节拍数， USpSp表示子矩阵①乘子矩阵④的总计算时钟节拍数，M表示矩阵A的行数，N表示矩 阵A的列数，R表示矩阵B的列数，S」表示矩阵B被按列划分的列数。
[0088] 情形 3 %= 1，k 2= 0):
[0089] 在这种情况下，与情形2类似，m= \jn=R/Sj，式（2)的评估函数可简化 为式（2-7)。
[0091] 式（2-7)中，TjSpSj表示子矩阵①乘子矩阵③的总计算时钟节拍数， T2,3(Si，Sp表示子矩阵②乘子矩阵③的总计算时钟节拍数，M表示矩阵A的行数，N表示矩 阵A的列数，R表示矩阵B的列数，5,表示矩阵A被按行划分的行数。
[0092] 情形 4(k1= l，k2= 1):
[0093] 令S'i=M-SiXm，S'」=R-SjXn，为了更好地简化讨论并且不失一般性，假设 S'i>S'j，式（2)的评估函数可简化为式（2-8)。
[0095] 式（2-8)中，USpSj)表示子矩阵①乘子矩阵③的总计算时钟节拍数， USi，sp表示子矩阵①乘子矩阵④的总计算时钟节拍数，T2,3(Si，sp表示子矩阵②乘子矩 阵③的总计算时钟节拍数，T2,4(Si，SP表示子矩阵②乘子矩阵④的总计算时钟节拍数。
[0096] 对于情形1至3,当Si的值小于或等于处理单元链中处理单元的个数时，评估函数 可取得最小值，因此在这三种情况下，最优的分块大小Si应当小于或等于处理单元链中处 理单元的个数P。对于情形4,无法直接求得评估函数的最小值，因此必须借助数学工具对 其进行分析，求得最优分块。对于具体的矩阵规模，M，N，R的值是已知的，也就意味着评估 函数是一元函数。本实施例中，具体是使用MATLAB数学工具对评估函数生成满足式（1)约 束的可视化图像，以矩阵A被按行划分的行数Si和矩阵B被按列划分的列数S/目等为前 提，可以发现评估函数的自变量（分块大小）只能取整数值，并且在整数点函数导数为零， 说明评估函数存在很多局部最值，本实施例通过作图分析的方式确定最优的分块大小 通过可视化图像确定评估函数定义域内的全局最小值，得到分块大小集合内最优的分块大 小Si〇
[0097] 以规模为M= 128、N= 1728、R= 169的矩阵乘为例，此时评估函数的表达式具体 如式（2-9)所示，并且自变量Si满足约束条件如式（2-10)所示：
[0099]
[0100] 式（2-9)中，51表示矩阵A被按行划分的行数，即最优的分块大小；L128/&」表示 对128/Si的结果向下取整；max为取较大值函数。本实施例中，具体是使用MATLAB数学工 具对评估函数生成满足式（1)约束的可视化图像如图4所示。从图4可知，最优的分块大 小Si在取值43时，评估函数达到最小值。子块大小Si的选取对矩阵加速器的性能发挥影 响很大，只有唯一的一个最优的分块大小使矩阵加速器发挥最高计算效率。
[0101] 本实施例中，步骤3)的详细步骤包括：
[0102] 3. 1)将矩阵A以规模为SjN的子块为单位按行分块，得到包含整数个规模为SjN 的子块的子矩阵①，如果仍有剩余的不规则的子块，则将剩余的不规则的子块作为子矩阵 ②；将矩阵B以规模为N*Sd^子块为单位按行分块，得到包含整数个规模为N*S,的子块的 子矩阵③，如果仍有剩余的不规则的子块，则将剩余的不规则的子块作为子矩阵④；将矩阵 A和矩阵B的矩阵乘运算视为得到的所有子矩阵的两两相乘，根据得到的子矩阵数量建立 使用指针相连接的DMA描述符链表，使得每一对相乘的子矩阵对应一个DMA描述符链表；
[0103] 3. 2)选择一个子块乘运算作为当前子块乘运算；
[0104] 3. 3)获取当前子块乘运算在矩阵A中对应子块数据的首地址、传输长度及传输步 长，获取当前子块乘运算在矩阵B中对应子块数据的首地址、传输长度及传输步长，将当前 子块乘运算在矩阵A中对应子块数据的首地址、传输长度及传输步长和当前子块乘运算在 矩阵B中对应子块数据的首地址、传输长度及传输步长封装生成一个DMA描述符，根据当前 子块乘运算所属的一对相乘的子矩阵确定对应的DMA描述符链表，并将该DMA描述符写入 对应的DMA描述符链表中；
[0105] 3. 4)判断是否所有子块乘运算已经完成处理，如果尚未完成所有子块乘运算的处 理，则选择下一个子块乘运算作为当前子块乘运算，跳转执行步骤3. 3);否则，如果已经完 成所有子块乘运算的处理则跳转执行步骤4)。
[0106] 假设矩阵A和矩阵B是两个规模为4*4的矩阵，当最优的子块大小Si= 2时，将矩 阵A分块为包含整数个规模为S#N的子块的子矩阵①、没有剩余的不规则的子矩阵②，将 矩阵B分块为包含整数个规模为N*Si的子块的子矩阵③、没有剩余的不规则的子矩阵④， 因此只存在一对相乘的子矩阵子矩阵①乘子矩阵③，此时矩阵共被划分为四个子块，矩阵 乘A*B总共需要计算四次分块乘法（A与B子块两两相乘）。每对子块乘法对应一个DMA描 述符（BufferDescriptor,BD)，因此一共有BD0~BD4共四个BD，四个BD组织成一条DMA 描述符链表存于主存中，如图5所示。DMA根据当前BD的信息就能够知道下一个BD在内 存的位置，因此DMA只需知道第一个BD的首地址即可完成对所有BD的读取，每个BD包含 了对应子块乘法的操作数的读取信息，DMA正是根据这些信息读取相应的数据。在整个读 数过程中，CPU无需参与，只需等待最终结果写回即可，这极大减少了主机和DMA的通信开 销。需要说明的是，如果出现子块数不是整数时，我们需要配置多个BD链，最差的情况下需 要配置四条BD链：最差的情况下，将矩阵A分块为包含整数个规模为SjN的子块的子矩阵 ①和剩余的不规则的子矩阵②、将矩阵B分块为包含整数个规模为N*Si的子块的子矩阵③ 和剩余的不规则的子矩阵④，因此共包含子矩阵①乘子矩阵③、子矩阵①乘子矩阵④、子矩 阵②乘子矩阵③、子矩阵②乘子矩阵④共四对相乘的子矩阵，此时则需要生成四个DMA描 述符链表，相邻DMA描述符链表之间使用指针相连接的DMA描述符链表。
[0107] 本实施例中，步骤4)的详细步骤包括：
[0108] 4. 1)CPU将DMA描述符链表中第一个DMA描述符的首地址配置给矩阵乘加速器的 DMA，矩阵乘加速器的DMA开始根据首地址读取第一个DMA描述符作为 当前描述符；
[0109] 4. 2)矩阵乘加速器的DMA解析当前描述符，得到当前描述符对应子块乘运算所需 的一对子块在外存中的地址及传输长度，并根据预设的传输步长和当前描述符中携带的一 对子块的地址及传输长度，以分时的方式交换读取当前描述符对应子块乘运算所需数据， 并将输入数据存入FIFO缓存；
[0110] 4.3)基于FIFO缓存中当前描述符对应子块乘运算所需数据，通过矩阵乘加速器 中至少一条处理单元链中的前Si个处理单元进行当前描述符对应子块乘运算；
[0111] 4. 4)将当前描述符对应子块乘运算的结果通过矩阵乘加速器的DMA分别写回主 存；
[0112] 4. 5)判断DMA描述符链表中的所有DMA描述符是否已经处理完毕，如果尚未处理 完毕，则从DMA描述符链表中选择下一个DMA描述符作为当前描述符，跳转执行步骤4. 2); 如果已经处理完毕，则结束并退出。
[0113] 本实施例中，步骤4. 3)的详细步骤包括：
[0114] 4. 3. 1)在矩阵乘加速器中确定一条用于当前描述符对应子块乘运算的处理单元 链，处理单元链至少包含Si个处理单元；矩阵乘加速器中设有用于控制各个处理单元工作 状态的状态机，状态机包括预取阶段、预取及计算阶段、提交阶段共三个状态；初始化变量 k为0,控制状态机进入预取状态；
[0115] 4. 3. 2)在预取阶段，针对当前描述符对应子块乘运算，矩阵乘加速器的DMA将属 于矩阵A的子块中的第k列数据附带上Sif处理单元的编号信息后发出给处理单元链，处 理单元链的前Si个处理单元根据编号信息预取编号对应的数据，控制状态机进入预取及计 算阶段；
[0116] 4. 3. 3)在预取及计算阶段，矩阵乘加速器的DMA将属于矩阵B的子块中第k行数 据以数据流的形式依次通过处理单元链的前Si个处理单元，且前Sif处理单元每一拍将接 收到的数据和属于矩阵A的子块中的第k列数据部分进行乘法运算，并将乘法运算结果和 片上存储中存储的上一次累加结果进行累加，将累加结果写回片上存储作为下一次累加的 输入；同时，处理单元链的前Sif处理单元分别根据编号信息预取属于矩阵A的子块中的 第k+1列数据；
[0117] 4. 3. 4)判断变量k的值是否等于N减1，如果变量k的值等于N减1，则控制状态 机进入提交阶段，跳转执行步骤4. 3. 5);否则，将初始化变量k加1，控制状态机进入预取及 计算阶段，跳转执行步骤4. 3. 3);
[0118] 4. 3. 5)在提交阶段下，处理单元链的前Si个处理单元分别将最后一次累加的结果 递给矩阵乘加速器的DMA，最终由矩阵乘加速器的DMA将结果写入主存中的指定区域。
[0119] 如图6所示，本实施例在步骤3)中，CPU根据矩阵信息生成各个子块乘运算对应 的DMA描述符（BufferDescriptor,BD)，BD组织成四条DMA描述符链表存于主存中；DMA 根据当前BD的信息就能够知道下一个BD在内存的位置，因此DMA只需知道第一个BD的首 地址即可完成对所有BD的读取，每个BD包含了对应子块乘法的操作数的读取信息，DMA正 是根据这些信息读取相应的数据。在整个读数过程中，CPU无需参与，只需等待最终结果写 回即可，这极大减少了主机和DMA的通信开销。在进入步骤4)后，首先CPU将BDO(第一条 DMA描述符链表的第一个DMA描述符）的首地址配置给DMA并启动DMA。DMA在工作状态 下，通过数据流的形式，自动读取DMA描述符链表中的DMA描述符至FIFO中，当读到第四个 DMA描述符链表BD4的最后一个DMA描述符时，DMA读取结束，等待矩阵加速器将结果写回 主存。
[0120] 如图7所示，本实施例中包含矩阵加速器的加速器系统基于xilinx公司的Zynq XC7Z045芯片上实现，主要包括处理机系统、片外存储和矩阵乘加速器，其中处理机系统包 括CPU、AXI总线模块和存储控制器，CPU具体采用ARM双核A9处理器，CPU可以通过AXI 总线模块经由存储控制器访问主存。处理机系统为矩阵加速器提供了AXI总线高速访存接 口，矩阵加速器通过该接口也能访问主存，并与CPU共享主存空间。矩阵加速器由DMA和处 理单元（ProcessingElement，PE)链（PE0~PEn)组成，每一个DMA和一个处理单元链相 连，且由链首的PE0直接和DMA交互，DMA直接与高速访存接口相连接，可直接从主存读取 加速器所需的数据。处理单元链是由若干个处理单元（PE)组成的链式结构，每个处理单元 的结构完全一致，数据在相邻的处理单元间进行传递。初始化的时候处理数据存储在片外 存储中，计算过程中的中间数据存在片上存储中，最终的结果又写回片外存储。具体过程是 第一个处理单元PE0从DMA读取数据，依次向后一个PE传递数据，每个PE保存中间结果， 等到全部计算完毕后，计算结果以相反方向向前一个PE传递，并最终由DMA将结果写回主 存。矩阵加速器可扩展为若干条处理单元链，每条处理单元链都有其对应的DMA。由于高性 能接口之间相对独立，都能够提供稳定的带宽，因此在DMA协助下，处理单元链之间可以高 效并行地完成计算任务。
[0121] DMA主要用以管理片外存储于加速器间的数据传输。它的主要特点是计算大规模 矩阵时往往只需要CPU配置一次，配置信息足够完整并存于片外存储中，DMA自动读取配置 信息后就能够完成所有的数据传输工作。参见图6,DMA的工作步骤主要包括：步骤1，CPU 根据矩阵分块后的矩阵的数据首地址、传输长度等信息生成DMA描述符，并将其组织成为 链表结构，存入片外存储。一个矩阵子块乘法对应一个描述符，描述符之间用指针相连接。 步骤2,CPU将描述符链表的在外存的首地址配置给DMA。DMA根据该地址自动读取DMA描 述符，然后解析出描述符的配置信息，得到矩阵输入数据在外存的地址，并根据传输步长， 传输长度等信息，以分时方式交换读取矩阵的输入数据。获得输入数据后，DMA将其存入相 应的FIFO缓存，为加速器提供操作数据。步骤3,DMA读取完一个子块乘法所需的操作数之 后，可以根据下个描述符的地址（当前描述符的信息）继续读取下一个子块乘法的操作数， 在此过程中，DMA还负责将上一个子块乘法的运算结果写回外存。步骤4,重复步骤3直到 DMA读取到最后一个DMA描述符。
[0122] 如图8所示，本实施例中每个处理单元由计算模块、计算数据传送模块、数据存储 模块和阶段同步控制模块组成。
[0123] 计算单元采用Xilinx公司提供的浮点运算IP核进行实现。计算单元包括乘法器 和加法器，其中加法器的流水线级别为11。
[0124] 计算数据传送模块由FIF0_A、FIF0_B、FIF0_C共三个FIFO队列组成，FIF0_A和 FIFO_B负责将输入数据传递到下一个PE。FIFO_C负责将结果传递到前一个PE中。
[0125] 数据存储模块主要包括一个双端口BRAM(MEM_C)和地址生成器。当乘法器产生第 一个结果时，地址生成器产生读信号和读地址，从MEM_C读取加法器的另一个操作数，也就 是上一阶段的临时数据送往加法器；当加法器的第一个结果产生时，地址生成器产生写信 号和写地址，结果被再次存入MEM_C。
[0126] 阶段同步控制模块内部主要实现了两个计数器，当矩阵A子块的新的一列和矩阵 B子块新的一行进入PE时，两个计数器分别开始计数。每个计数器都有相应的阈值，一旦 其中一个计数器达到阈值，阶段同步控制模块就会使当前PE停止从前一个PE读取相对应 的数据。这样一来达到阈值的计数器就会停止计数，一直到另一个计数器也达到其阈值， 两个计数器同时被清零，重新开始为下一阶段的数据计数。通过阶段同步保证了矩阵A的 某一列的第一个数据总是和矩阵B对应行的第一个数据同时进入每个PE。每个PE都有一 个ID(PID)，矩阵A的每个数据在进入矩阵链之前被加入了编号信息用以表征这个数据归 属于哪个处理单元。阶段同步控制模块，用以解决支持可变分块引起的数据不同步问题；此 外，每个处理单元还包括用以控制与相邻处理单元间的数据交互的控制逻辑。
[0127] 需要说明的是，矩阵加速器中处理单元的结构实现并不局限于如图8所示的特定 结构，毫无疑问，本领域技术人员也可以根据需要采用其他结构的处理单元来实现矩阵加 速器中的子块乘法运算，故在此不再赘述。
[0128] 根据图2的原理可知，编号信息依次为0，1，…，Sri.iS'P时，那些满足PID> =PE不允许启动工作。我们通过对数据编号和PID进行对比，使得当数据编号小于 当前PID时，数据才允许被写入FIFO_A并传递到下一个PE。这样一来，数据就只传递到前 SifPE中，我们的矩阵乘法器就能支持可变大小的分块。
[0129] 本实施例中，每个处理单元内部实现了 一个状态机用来控制数据流。处理单元工 作时主要有三个阶段，我们以处理单元链中的第一个处理单元（PE0)为例：第一，预取阶 段。PE0从DMA的FIFO中读取矩阵A的第一列数据，并通过FIFO_A传递给相邻的PE。数 据从FIF0_A出来后在高位加入了编号信息，每个PE根据自己的ID号预取各自的数据。第 二，计算-预取阶段。PE0开始同时读取矩阵A的第k+1列以及矩阵B的第k行数据并通过 FIF0_A和FIF0_B传递给相邻的PE;矩阵B的数据与预取的A数据相乘，同时更新A的另一 个预取缓冲（双缓冲，TA0，TA1)。经过一定延迟后对乘法结果进行加法运算，其中另一个操 作数来自MEM_C，加法计算完毕后，又将结果存入MEM_C。第三阶段，提交阶段。每个PE将 最后一次累加结果存入FIF0_C，同时MEM_C内容清零。存入完毕后 ，每个PE从后一个PE读 取结果存入FIF0_C，由此结果数据依次前移至DMA的FIF0_C中。
[0130] 本实施例以典型卷积神经网络（CNN)中如表1所示的五个非均匀矩阵规模 Convl~Conv5为例，应用本实施例以及经典分块算法分别对表1中的五个规模的矩阵乘进 行加速对比，表2是本实施例中每个卷积层的矩阵规模通过优化分块策略算出的最优的分 块大小Sp
[0131] 表1 :五个非均匀矩阵规模实例表。
[0133] 表2 :五个非均匀矩阵规模实例的最优的分块大小31数据表。
[0135] 最终，得到的性能对比图如图9所示。参见图9可知，当矩阵规模不均匀程度比较 高时（矩阵的行列大小相差较大），本实施例的矩阵乘甲酸方法能达到很高的性能提升，例 如对于Conv5,本实施例的矩阵乘甲酸方法相对于经典分块算法有大约12%的性能提升。
[0136] 以上所述仅是本发明的优选实施方式，本发明的保护范围并不仅局限于上述实施 例，凡属于本发明思路下的技术方案均属于本发明的保护范围。应当指出，对于本技术领域 的普通技术人员来说，在不脱离本发明原理前提下的若干改进和润饰，这些改进和润饰也 应视为本发明的保护范围。
【主权项】
1. 一种支持可变分块的矩阵乘加速方法，其特征在于步骤包括： 1) 输入矩阵乘运算所需的M*N的矩阵A和N*R的矩阵B ; 2) 根据矩阵A和矩阵B的规模确定子块大小Si,将矩阵A以规模为SjN的子块为单位 进行按行分块，将矩阵B以规模为N*Si的子块为单位进行按列分块，使得矩阵乘运算等同 为多个子块乘运算； 3) 为每一个子块乘运算所需数据生成一个DM描述符，将所有子块乘运算的DM描述 符构建DM描述符链表并存入主存； 4) 针对每一个子块乘运算，通过矩阵乘加速器的DM从主存的DM描述符链表读取子 块乘运算所需数据，然后通过矩阵乘加速器中至少一条处理单元链中的前Si个处理单元进 行子块乘运算，并将各个子块乘运算的结果通过矩阵乘加速器的DM分别写回主存。2. 根据权利要求1所述的支持可变分块的矩阵乘加速方法，其特征在于，所述步骤2) 中根据矩阵A和矩阵B的规模确定子块大小详细步骤包括： 2. 1)输入矩阵A和矩阵B的矩阵规模，所述矩阵规模包括矩阵A的行数M、矩阵A的列 数N、矩阵B的列数R三者的值； 2. 2)根据矩阵规模获取满足式（1)所示约束的分块大小集合；式⑴中，示矩阵A被按行划分的行数，S廣示矩阵B被按列划分的列数，M表 示矩阵A的行数，N表示矩阵A的列数，BW表示矩阵乘加速器中单条处理单元链的访 存带宽，HiaxiSi, SjI表示从SjP S」中取较大值，F表示访存频率；P表示矩阵乘加速器 中单条处理单元链的处理单元数量，Stageadd表示矩阵乘加速器中的加法器流水级数， Hiaxm-SiXm, R-SjXn}表示从M-SiXm和R-SjXn中取较大值，R表示矩阵B的列数，m表示 矩阵A的行数M除以矩阵A被按行划分的行数结果的向下取整值，η表示矩阵B的列 数R除以矩阵B被按列划分的列数\_的结果的向下取整值； 2. 3)将矩阵A视为包含整数个规模为S#N的子块的子矩阵①和剩余的不规则的子矩 阵②、将矩阵B视为包含整数个规模为N*Sj^子块的子矩阵③和剩余的不规则的子矩阵 ④，将矩阵A和矩阵B的矩阵乘运算视为子矩阵①~子矩阵④四者中的两两相乘，建立式 (2)所示的评估函数； f (Si, Sj) = 8,.)+^3(8^ Sj) (2) 式（2)中，f (Si, Sp表示评估函数的值，51表示矩阵A被按行划分的行数，S ^表示矩阵 B被按列划分的列数，表示子矩阵①乘子矩阵③的计算时钟节拍数，Tli4(SilSj) 表示子矩阵①乘子矩阵④的计算时钟节拍数，T2i3(Si，sp表示子矩阵②乘子矩阵③的计算 时钟节拍数，(Si, sp表示子矩阵②乘子矩阵④的计算时钟节拍数； 2.4)使用MATLAB数学工具对式（2)所示的评估函数生成满足式（1)约束的可视化图 像，以矩阵A被按行划分的行数Si和矩阵B被按列划分的列数S d目等为前提，通过所述可 视化图像确定评估函数定义域内的全局最小值，得到分块大小集合内最优的分块大小Sp3.根据权利要求2所述的支持可变分块的矩阵乘加速方法，其特征在于，所述步骤3) 的详细步骤包括： 3. 1)将矩阵A以规模为S#N的子块为单位按行分块，得到包含整数个规模为S#N的 子块的子矩阵①，如果仍有剩余的不规则的子块，则将剩余的不规则的子块作为子矩阵②； 将矩阵B以规模为N*Sj^子块为单位按行分块，得到包含整数个规模为N*S i的子块的子矩 阵③，如果仍有剩余的不规则的子块，则将剩余的不规则的子块作为子矩阵④；将矩阵A和 矩阵B的矩阵乘运算视为得到的所有子矩阵的两两相乘，根据得到的子矩阵数量建立使用 指针相连接的DM描述符链表，使得每一对相乘的子矩阵对应一个DM描述符链表； 3. 2)选择一个子块乘运算作为当前子块乘运算； 3. 3)获取当前子块乘运算在矩阵A中对应子块数据的首地址、传输长度及传输步长， 获取当前子块乘运算在矩阵B中对应子块数据的首地址、传输长度及传输步长，将当前子 块乘运算在矩阵A中对应子块数据的首地址、传输长度及传输步长和当前子块乘运算在矩 阵B中对应子块数据的首地址、传输长度及传输步长封装生成一个DM描述符，根据当前子 块乘运算所属的一对相乘的子矩阵确定对应的DM描述符链表，并将该DM描述符写入对 应的DM描述符链表中； 3. 4)判断是否所有子块乘运算已经完成处理，如果尚未完成所有子块乘运算的处理， 则选择下一个子块乘运算作为当前子块乘运算，跳转执行步骤3. 3);否则，如果已经完成 所有子块乘运算的处理则跳转执行步骤4)。4. 根据权利要求3所述的支持可变分块的矩阵乘加速方法，其特征在于，所述步骤4) 的详细步骤包括： 4. I) CPU将DM描述符链表中第一个DM描述符的首地址配置给矩阵乘加速器的DMA， 矩阵乘加速器的DM开始根据首地址读取第一个DM描述符作为当前描述符； 4. 2)矩阵乘加速器的DM解析当前描述符，得到当前描述符对应子块乘运算所需的一 对子块在外存中的地址及传输长度，并根据预设的传输步长和当前描述符中携带的一对子 块的地址及传输长度，以分时的方式交换读取当前描述符对应子块乘运算所需数据，并将 输入数据存入FIFO缓存； 4. 3)基于FIFO缓存中当前描述符对应子块乘运算所需数据，通过矩阵乘加速器中至 少一条处理单元链中的前Si个处理单元进行当前描述符对应子块乘运算； 4. 4)将当前描述符对应子块乘运算的结果通过矩阵乘加速器的DM分别写回主存； 4. 5)判断DMA描述符链表中的所有DMA描述符是否已经处理完毕，如果尚未处理完毕， 则从DM描述符链表中选择下一个DM描述符作为当前描述符，跳转执行步骤4. 2);如果 已经处理完毕，则结束并退出。5. 根据权利要求4所述的支持可变分块的矩阵乘加速方法，其特征在于，所述步骤 4. 3)的详细步骤包括： 4. 3. 1)在矩阵乘加速器中确定一条用于当前描述符对应子块乘运算的处理单元链，所 述处理单元链至少包含Si个处理单元；所述矩阵乘加速器中设有用于控制各个处理单元工 作状态的状态机，所述状态机包括预取阶段、预取及计算阶段、提交阶段共三个状态；初始 化变量k为0,控制状态机进入预取状态； 4. 3. 2)在预取阶段，针对当前描述符对应子块乘运算，矩阵乘加速器的DM将属于矩 阵A的子块中的第k列数据附带上SiA处理单元的编号信息后发出给处理单元链，处理单 元链的前Si个处理单元根据编号信息预取编号对应的数据，控制状态机进入预取及计算阶 段； 4. 3. 3)在预取及计算阶段，矩阵乘加速器的DM将属于矩阵B的子块中第k行数据以 数据流的形式依次通过处理单元链的前Si个处理单元，且前S i个处理单元每一拍将接收 到的数据和属于矩阵A的子块中的第k列数据部分进行乘法运算，并将乘法运算结果和片 上存储中存储的上一次累加结果进行累加，将累加结果写回片上存储作为下一次累加的输 入；同时，处理单元链的前SiA处理单元分别根据编号信息预取属于矩阵A的子块中的第 k+Ι列数据； 4. 3. 4)判断变量k的值是否等于N减1，如果变量k的值等于N减1，则控制状态机进 入提交阶段，跳转执行步骤4. 3. 5);否则，将初始化变量k加1，控制状态机进入预取及计算 阶段，跳转执行步骤4.3.3); 4. 3. 5)在提交阶段下，处理单元链的前SiA处理单元分别将最后一次累加的结果递给 矩阵乘加速器的DMA，最终由矩阵乘加速器的DMA将结果写入主存中的指定区域。
【专利摘要】本发明公开了一种支持可变分块的矩阵乘加速方法，步骤包括：输入矩阵A和矩阵B，根据矩阵A和矩阵B的规模确定子块大小Si，将矩阵A以规模为Si*N的子块为单位进行按行分块，将矩阵B以规模为N*Si的子块为单位进行按列分块，为每一个子块乘运算所需数据生成一个DMA描述符，将所有DMA描述符构建DMA描述符链表；针对每一个子块乘运算，根据主存的DMA描述符链表读取子块乘运算所需数据，通过矩阵乘加速器中的一条处理单元链进行子块乘运算，并将结果通过DMA写回主存。本发明具有支持可变分块，可根据分块大小调整使用的处理单元数目，加速非均匀矩阵乘运算加速效率高的优点。
【IPC分类】G06F7/523, G06F17/16
【公开号】CN104899182
【申请号】CN201510312188
【发明人】文梅, 沈俊忠, 乔寓然, 杨乾明, 苏华友, 肖涛, 陈照云, 张春元
【申请人】中国人民解放军国防科学技术大学
【公开日】2015年9月9日
【申请日】2015年6月9日
【技术领域】
[0001] 本发明涉及嵌入式平台下的矩阵乘加速技术，具体涉及一种支持可变分块的矩阵 乘加速方法。
【背景技术】
[0002] 随着半导体制造工艺的发展以及集成电路技术的进步，单芯片上能够集成越来越 多的晶体管，使用可编程器件，特别是FPGA(FieldProgrammableGateArray)芯片进行设 计成为了目前构建嵌入式系统以及硬件加速平台的一种重要方式。当前FPGA芯片提供了 专用的算术模块、大量的逻辑资源和存储资源，以及外部存储器接口、网络接口和其它外围 接口，为构建高性能计算系统提供了条件，也使FPGA可重构计算系统成为加速科学计算的 一种重要选择。当前的FPGA芯片相对于DSP(DigitalSignalProcessor,数字信号处理 器）具有可编程的优势，同时能够并行处理海量数据，既具有通用处理器的灵活性，又具有 ASIC(ApplicationSpecificIntegratedCircuit,专用集成电路）的高性能，在嵌入式计 算领域备受青睐。
[0003] 浮点矩阵乘法是数字信号的基本算法，同时也是许多科学计算方法的基本运算。 在数字图像处理，计算机视觉的快速处理以及工业实时控制等领域都被广泛应用。但由于 实际应用中，浮点矩阵规模通常较大，矩阵乘算法本身复杂度较高、处理效率较低，成为限 制系统性能提升的瓶颈所在，因此为此类应用设计高性能的硬件结构是当前FPGA结构设 计的研宄热点。
[0004] 尽管近年来提出了很多基于FPGA的矩阵乘加速器设计，但是都缺乏对非均匀的 大规模矩阵加速的讨论和支持，这种大规模矩阵的特征是行列数相差很大（> =10倍），并 且广泛存在于很多现代应用领域中，如图像处理，深度学习等。在这些应用领域之中，矩阵 乘占据了计算量的主要部分。由于单片FPGA芯片上存储资源和计算资源十分有限，加速大 规模矩阵乘时往往需要对矩阵进行分块。对于链式结构的矩阵乘加速器，虽然对大多数大 规模矩阵的加速效果十分明显，但对于加速非均匀矩阵时的计算效率却很低，其主要原因 就是这类加速器往往只支持固定的分块，也就是说，分块大小与矩阵链长（矩阵链中处理 单元个数）相等或者是其倍数。当分块大小与矩阵加速器链长不匹配时，加速器计算效率 会发生明显的下降。据我们所知，至今还没有公开文献涉及支持可变分块的矩阵乘法器设 计，也没有经典分块算法基础上的分块优化问题方面的相关研宄，因此如何选择最优的分 块，使得矩阵乘加速器获得更高的计算效率，以更好地适应现代应用的加速需求，已经成为 一项亟待解决的关键技术问题。
【发明内容】
[0005] 本发明要解决的技术问题是：针对现有技术的上述技术问题，提供一种支持可变 分块，可根据分块大小调整使用的处理单元数目，加速非均匀矩阵乘运算加速效率高的支 持可变分块的矩阵乘加速方法。
[0006] 为了解决上述技术问题，本发明采用的技术方案为：
[0007] -种支持可变分块的矩阵乘加速方法，步骤包括：
[0008] 1)输入矩阵乘运算所需的M*N的矩阵A和N*R的矩阵B;
[0009] 2)根据矩阵A和矩阵B的规模确定子块大小Si，将矩阵A以规模为Si*N的子块为 单位进行按行分块，将矩阵B以规模为N*Si的子块为单位进行按列分块，使得矩阵乘运算 等同为多个子块乘运算；
[0010] 3)为每一个子块乘运算所需数据生成一个DMA描述符，将所有子块乘运算的DMA描述符构建DMA描述符链表并存入主存；
[0011] 4)针对每一个子块乘运算，通过矩阵乘加速器的DMA从主存的DMA描述符链表读 取子块乘运算所需数据，然后通过矩阵乘加速器中至少一条处理单元链中的前Sif处理单 元进行子块乘运算，并将各个子块乘运算的结果通过矩阵乘加速器的DMA分别写回主存。
[0012] 优选地，所述步骤2)中根据矩阵A和矩阵B的规模确定子块大小Si的详细步骤 包括：
[0013] 2. 1)输入矩阵A和矩阵B的矩阵规模，所述矩阵规模包括矩阵A的行数M、矩阵A 的列数N、矩阵B的列数R三者的值；
[0014] 2. 2)根据矩阵规模获取满足式⑴所示约束的分块大小集合；
[0016] 式（1)中，3,表示矩阵A被按行划分的行数，S」表示矩阵B被按列划分的列数， M表示矩阵A的行数，N表示矩阵A的列数，BW表示矩阵乘加速器中单条处理单元链的 访存带宽，max以，Sj表示从SJP\中取较大值，F表示访存频率；P表示矩阵乘加速 器中单条处理单元链的处理单元数量，Stageadd表示矩阵乘加速器中的加法器流水级数， maxm-SiXi^R-SjXn}表示从M-SiXm和R-SjXn中取较大值，R表示矩阵B的列数，m表示 矩阵A的行数M除以矩阵A被按行划分的行数结果的向下取整值，n表示矩阵B的列 数R除以矩阵B被按列划分的列数&的结果的向下取整值；
[0017] 2. 3)将矩阵A视为包含整数个规模为S#N的子块的子矩阵①和剩余的不规则的 子矩阵②、将矩阵B视为包含整数个规模为N*Si的子块的子矩阵③和剩余的不规则的子矩 阵④，将矩阵A和矩阵B的矩阵乘运算视为子矩阵①~子矩阵④四者中的两两相乘，建立式 (2)所示的评估函数；
[0018]f(S^ Sj)=T1； 3(S^ Sj)+T1； 4(S^ Sj)+T2； 3(S^ Sj)+T2； 4(S^ Sj)(2)
[0019] 式⑵中，f(Si，Sp表示评估函数的值，5,表示矩阵A被按行划分的行数，5」表 示矩阵B被按列划分的列数，USpSp表示子矩阵①乘子矩阵③的计算时钟节拍数， USi，Sj表示子矩阵①乘子矩阵④的计算时钟节拍数，T2,3(Si，Sj表示子矩阵②乘子矩阵 ③的计算时钟节拍数，T2,4(Si，SP表示子矩阵②乘子矩阵④的计算时钟节拍数；
[0020] 2. 4)使用MATLAB数学工具对式（2)所示的评估函数生成满足式（1)约束的可视 化图像，以矩阵A被按行划分的行数Si和矩阵B被按列划分的列数Sd目等为前提，通过所 述可视化图像确定评估函数定义域内的全局最小值，得到分块大小集合内最优的分块大小 si〇
[0021] 优选地，所述步骤3)的详细步骤包括：
[0022] 3. 1)将矩阵A以规模为SjN的子块为单位按行分块，得到包含整数个规模为SjN 的子块的子矩阵①，如果仍有剩余的不规则的子块，则将剩余的不规则的子块作为子矩阵 ②；将矩阵B以规模为N*Sd^子块为单位按行分块，得到包含整数个规模为N*S,的子块的 子矩阵③，如果仍有剩余的不规则的子块，则将剩余的不规则的子块作为子矩阵④；将矩阵 A和矩阵B的矩阵乘运算视为得到的所有子矩阵的两两相乘，根据得到的子矩阵数量建立 使用指针相连接的DMA描述符链表，使得每一对相乘的子矩阵对应一个DMA描述符链表；
[0023] 3. 2)选择一个子块乘运算作为当前子块乘运算；
[0024] 3. 3)获取当前子块乘运算在矩阵A中对应子块数据的首地址、传输长度及传输步 长，获取当前子块乘运算在矩阵B中对应子块数据的首地址、传输长度及传输步长，将当前 子块乘运算在矩阵A中对应子块数据的首地址、传输长度及传输步长和当前子块乘运算在 矩阵B中对应子块数据的首地址、传输长度及传输步长封装生成一个DMA描述符，根据当前 子块乘运算所属的一对相乘的子矩阵确定对应的DMA描述符链表，并将该DMA描述符写入 对应的DMA描述符链表中；
[0025] 3. 4)判断是否所有子块乘运算已经完成处理，如果尚未完成所有子块乘运算的处 理，则选择下一个子块乘运算作为当前子块乘运算，跳转执行步骤3. 3);否则，如果已经完 成所有子块乘运算的处理则跳转执行步骤4)。
[0026] 优选地，所述步骤4)的详细步骤包括：
[0027] 4. 1)CPU将DMA描述符链表中第一个DMA描述符的首地址配置给矩阵乘加速器的 DMA，矩阵乘加速器的DMA开始根据首地址读取第一个DMA描述符作为当前描述符；
[0028] 4. 2)矩阵乘加速器的DMA解析当前描述符，得到当前描述符对应子块乘运算所需 的一对子块在外存中的地址及传输长度，并根据预设的传输步长和当前描述符中携带的一 对子块的地址及传输长度，以分时的方式交换读取当前描述符对应子块乘运算所需数据， 并将输入数据存入FIFO缓存；
[0029] 4. 3)基于FIFO缓存中当前描述符对应子块乘运算所需数据，通过矩阵乘加速器 中至少一条处理单元链中的前Si个处理单元进行当前描述符对应子块乘运算；
[0030] 4. 4)将当前描述符对应子块乘运算的结果通过矩阵乘加速器的DMA分别写回主 存；
[0031] 4. 5)判断DMA描述符链表中的所有DMA描述符是否已经处理完毕，如果尚未处理 完毕，则从DMA描述符链表中选择下一个DMA描述符作为当前描述符，跳转执行步骤4. 2); 如果已经处理完毕，则结束并退出。
[0032] 优选地，所述步骤4. 3)的详细步骤包括：
[0033] 4. 3. 1)在矩阵乘加速器中确定一条用于当前描述符对应子块乘运算的处理单元 链，所述处理单元链至少包含Si个处理单元；所述矩阵乘加速器中设有用于控制各个处理 单元工作状态的状态机，所述状态机包括预取阶段、预取及计算阶段、提交阶段共三个状 态；初始化变量k为0,控制状态机进入预取状态；
[0034] 4. 3. 2)在预取阶段，针对当前描述符对应子块乘运算，矩阵乘加速器的DMA将属 于矩阵A的子块中的第k列数据附带上Sif处理单元的编号信息 后发出给处理单元链，处 理单元链的前Si个处理单元根据编号信息预取编号对应的数据，控制状态机进入预取及计 算阶段；
[0035] 4. 3. 3)在预取及计算阶段，矩阵乘加速器的DMA将属于矩阵B的子块中第k行数 据以数据流的形式依次通过处理单元链的前Si个处理单元，且前Sif处理单元每一拍将接 收到的数据和属于矩阵A的子块中的第k列数据部分进行乘法运算，并将乘法运算结果和 片上存储中存储的上一次累加结果进行累加，将累加结果写回片上存储作为下一次累加的 输入；同时，处理单元链的前Sif处理单元分别根据编号信息预取属于矩阵A的子块中的 第k+1列数据；
[0036] 4. 3. 4)判断变量k的值是否等于N减1，如果变量k的值等于N减1，则控制状态 机进入提交阶段，跳转执行步骤4. 3. 5);否则，将初始化变量k加1，控制状态机进入预取及 计算阶段，跳转执行步骤4. 3. 3);
[0037] 4. 3. 5)在提交阶段下，处理单元链的前Si个处理单元分别将最后一次累加的结果 递给矩阵乘加速器的DMA，最终由矩阵乘加速器的DMA将结果写入主存中的指定区域。
[0038] 本发明支持可变分块的矩阵乘加速方法具有下述优点：本发明通过确定子块大小 Sp将矩阵A以规模为S#N的子块为单位进行按行分块，将矩阵B以规模为N*Si的子块为单 位进行按列分块，使得矩阵乘运算等同为多个子块乘运算，为每一个子块乘运算所需数据 生成一个DMA描述符，将所有子块乘运算的DMA描述符构建DMA描述符链表并存入主存，针 对每一个子块乘运算，通过矩阵乘加速器的DMA从主存的DMA描述符链表读取子块乘运算 所需数据，然后通过矩阵乘加速器中至少一条处理单元链中的前Sif处理单元进行子块乘 运算，并将各个子块乘运算的结果通过矩阵乘加速器的DMA分别写回主存，增加了对可变 分块的支持，能够根据具体的分块大小调整使用的处理单元数目，通过确定子块大小3 1来 基于经典分块算法上的优化分块策略，当加速非均匀矩阵时，通过计算并应用最优分块能 够使得矩阵加速器具有很高的计算效率，具有支持可变分块，可根据分块大小调整使用的 处理单元数目，加速非均匀矩阵乘运算加速效率高的优点。
【附图说明】
[0039] 图1为本发明实施例的基本流程示意图。
[0040] 图2为本发明实施例将进行矩阵乘运算等同为多个子块乘运算的原理示意图。 [0041] 图3为本发明实施例将进行矩阵乘运算划分为子矩阵相乘的原理示意图。
[0042] 图4为本发明实施例中使用MATLAB数学工具生成的评估函数可视化图像。
[0043] 图5为本发明实施例中生成的一个DMA描述符链表的结构不意图。
[0044] 图6为本发明实施例中矩阵加速器的DMA的工作流程示意图。
[0045] 图7为本发明实施例中加速器系统的框架结构示意图。
[0046] 图8为本发明实施例中处理单元的数据流向结构示意图。
[0047] 图9为本发明实施例和采用经典分块算法的性能对比图。
【具体实施方式】
[0048] 如图1所示，本实施例支持可变分块的矩阵乘加速方法的步骤包括：
[0049] 1)输入矩阵乘运算所需的M*N的矩阵A和N*R的矩阵B;
[0050] 2)根据矩阵A和矩阵B的规模确定子块大小Si，将矩阵A以规模为Si*N的子块为 单位进行按行分块，将矩阵B以规模为N*Si的子块为单位进行按列分块，使得矩阵乘运算 等同为多个子块乘运算；
[0051] 3)为每一个子块乘运算所需数据生成一个DMA描述符，将所有子块乘运算的DMA 描述符构建DMA描述符链表并存入主存；
[0052] 4)针对每一个子块乘运算，通过矩阵乘加速器的DMA从主存的DMA描述符链表读 取子块乘运算所需数据，然后通过矩阵乘加速器中至少一条处理单元链中的前Sif处理单 元进行子块乘运算，并将各个子块乘运算的结果通过矩阵乘加速器的DMA分别写回主存。
[0053] 如图2所示，对于M*N的矩阵A和N*R的矩阵B的矩阵乘运算而言，其矩阵乘运算 得到的矩阵C的规模为M*R。因此，将矩阵A和N*R的矩阵B的矩阵乘运算划分为多个小 矩阵相乘，其中矩阵A按行划分为多个S#N子块，矩阵B按列划分为多个N*Sj子块，每对小 子块相乘就得到一个SjSj勺矩阵，也就是最后的结果C矩阵的一个子块。参见图2,现考 虑子块相乘的细节，可把A的S#N子块看成是由N个维度为列向量组成的矩阵，B的 N*Sj子块是由N个维度为S」的行向量组成的矩阵，记第k个列向量和行向量分别为Vk和 Uk(k〈 =N)，把Uk和Vk相乘的结果进行累加（N次）就得到矩阵C的SjSj子块。假设矩阵 乘加速器有Si个处理单元，每个处理单元每次缓冲Vk中的一个操作数，Uk中的所有操作数 依次进入每个处理单元，处理单元以流水线的形式高效地进行乘加运算，并且将加法结果 存入临时存储区（下一次加法运算要从存储区读取这一次的结果进行累加），这样一来，每 个处理单元并行地计算C的某一行的结果，这就是本实施例算法的核心思想。事实上当Si =SJt，矩阵乘加速器进行阶段同步所导致的空拍时最少，这种情况下矩阵乘加速器的性 能应该是最优的，因此本实施例仅考虑Si=Sj勺情况。
[0054] 本实施例中，步骤2)中根据矩阵A和矩阵B的规模确定子块大小Si的详细步骤 包括：
[0055] 2. 1)输入矩阵A和矩阵B的矩阵规模，矩阵规模包括矩阵A的行数M、矩阵A的列 数N、矩阵B的列数R三者的值；
[0056] 2. 2)根据矩阵规模获取满足式⑴所示约束的分块大小集合；
[0058] 式（1)中，31表示矩阵A被按行划分的行数，5』表示矩阵B被按列划分的列数， M表示矩阵A的行数，N表示矩阵A的列数，BW表示矩阵乘加速器中单条处理单元链的 访存带宽，max以，Sj表示从SJP\中取较大值，F表示访存频率；P表示矩阵乘加速 器中单条处理单元链的处理单元数量，Stageadd表示矩阵乘加速器中的加法器流水级数， maxm-SiXi^R-SjXn}表示从M-SiXm和R-SjXn中取较大值，R表示矩阵B的列数，m表示 矩阵A的行数M除以矩阵A被按行划分的行数结果的向下取整值，n表示矩阵B的列 数R除以矩阵B被按列划分的列数&的结果的向下取整值；
[0059] 2. 3)将矩阵A视为包含整数个规模为S#N的子块的子矩阵①和剩余的不规则的 子矩阵②、将矩阵B视为包含整数个规模为N*Si的子块的子矩阵③和剩余的不规则的子矩 阵④，将矩阵A和矩阵B的矩阵乘运算视为子矩阵①~子矩阵④四者中的两两相乘，建立式 (2)所示的评估函数；
[0060] f(S^Sj) =T1； 3 (S^Sj) +T1； 4 (S^Sj) +T2； 3 (S^Sj) +T2； 4 (S^Sj) (2)
[0061] 式（2)中，表示评估函数的值，51表示矩阵A被按行划分的行数，Sj-表 示矩阵B被按列划分的列数，USpSp表示子矩阵①乘子矩阵③的计算时钟节拍数， USi，Sj表示子矩阵①乘子矩阵④的计算时钟节拍数，T2,3(Si，Sj表示子矩阵②乘子矩阵 ③ 的计算时钟节拍数，T2,4(Si，SP表示子矩阵②乘子矩阵④的计算时钟节拍数；
[0062] 2. 4)使用MATLAB数学工具对式（2)所示的评估函数生成满足式（1)约束的可视 化图像，以矩阵A被按行划分的行数Si和矩阵B被按列划分的列数Sd目等为前提，通过可 视化图像确定评估函数定义域内的全局最小值，得到分块大小集合内最优的分块大小Sp
[0063] 如图3所示，将矩阵A和矩阵B的矩阵乘运算视为子矩阵①~子矩阵④四者中的 两两相乘后，子矩阵①包含整数个规模为SjN的子块，子矩阵③包含整数个规模为N*SJ9 子块，剩余的不规则的子矩阵②包含的为不规则额的子块（行数小于SJ，剩余的不规则的 子矩阵④包含的为不规则的子块（列数小于SJ，子矩阵①~子矩阵④四者中的两两相乘 即为子矩阵①乘子矩阵③、子矩阵①乘子矩阵④、子矩阵②乘子矩阵③、子矩阵②乘子矩阵 ④ 。本实施例中，由于子块乘运算要求在开始的时候预取Si个矩阵A的数据（子块第一 列），这个过程需要Si个节拍，每个PE计算SJSj需要花费maxtSpSjlXN个节拍。由于 矩阵乘加速器中的处理单元组织成链式结构，因此最后一个计算的处理单元相对于第一个 处理单元有Si个节拍的延迟，即整个处理单元流水线的建立时间。计算S，5」子块的结 果需要从主存读取（SJSPXN个数据，该过程需要时间（Si+SpXN/BWJW是单条链的访存 带宽。根据所有处理单元总计算时间大于取数时间，有以下不等式（1-1)成立。
[0064] (Sj+Sj)XN/BW^ (2XSj+max{S^Sj}XN)/F(1-1)
[0065] 式（1-1)中，3,表示矩阵A被按行划分的行数，S」表示矩阵B被按列划分的列数， N表示矩阵A的列数，BW表示矩阵乘加速器中单条处理单元链的访存带宽，max以，Sj表示 从SJP\中取较大值，F表示访存频率。
[0066] 此外，矩阵A被按行划分的行数Sp矩阵B被按列划分的列数&还有以下约束：首 先，根据算法，矩阵A被按行划分的行数Si的值不能大于矩阵乘加速器中一条处理单元链 中处理单元的数目；其次，考虑到流水线中的数据冲突，也就是当加法结果尚未写入片上存 储器时，加法器又需要从存储器读取该结果时，加法器就会读到旧值。因此，必须同时满足 式（1-2)所示的几个约束条件。
[0068] 式（1-2)中，51表示矩阵A被按行划分的行数，P表示矩阵乘加速器中单条处理单 元链的 处理单元数量，max以，S」}表示从SjPS」中取较大值，Stageadd表示矩阵乘加速器 中的加法器流水级数，11^復-5'111，1^，11}表示从1-5'111和1^，11中取较大值。本实 施例中，Stageadd表示矩阵乘加速器中的加法器流水级数具体值为11。结合前述的式（1-1) 和式（1-2)，即可推导得出本实施例前述的式（1)所示约束的分块大小集合。
[0069] 假设?q」，" =[RJ(即m表示矩阵A的行数M除以矩阵A被按行划分的行 数结果的向下取整值，n表示矩阵B的列数R除以矩阵B被按列划分的列数\的结果 的向下取整值），则有子矩阵①乘子矩阵③的总计算时钟节拍数USpSj如式（2-1)所 示，子矩阵①乘子矩阵④的总计算时钟节拍数USpSj如式（2-2)所示，子矩阵②乘子矩 阵③的总计算时钟节拍数T2,3(Si，Sp如式（2-3)所示，子矩阵②乘子矩阵④的总计算时钟 节拍数如式（2-4)所示，
[0070] 1\，3以，Sj) =mXnX(Si+maxtSi，SjXN) (2-1)
[0071] Sj) =k2XmX(maxtSj,R-SjXn}XN) (2-2)
[0072] T2； 3 (S^Sj) =k:XnX(max{M-SjXm,Sj}XN) (2-3)
[0073] T2； 4 (S^Sj) =k:Xk2X(max{M-SjXm,R-SjXn}XN) (2-4)
[0074] 式（2-1)~（2-4)中，maxtSi，S」}表示从SJPS」中取较大值，max{Si，R-SjXn} 表示从SJPR-S」Xn中取较大值，max{M-SiXm，S」}表示从M-SiXm和Sj中取较大值， maxm-SiXi^R-SjXn}表示从M-SiXm和R-SjXn中取较大值，m表示矩阵A的行数M除以 矩阵A被按行划分的行数结果的向下取整值，n表示矩阵B的列数R除以矩阵B被按 列划分的列数h的结果的向下取整值，S,表示矩阵A被按行划分的行数，S^表示矩阵B被 按列划分的列数，M表示矩阵A的行数，N表示矩阵A的列数，R表示矩阵B的列数。
[0075] 式（2-2)~（2-4)中，变量&和1^2的值由式（1-2)和（1-3)决定。
[0078] 式（1-2)和（1-3)中各个参数的定义可参见式（2-1)~（2-4)，在此不再赘述。
[0079] 参见式（2-1)~（2-4)可知，评估函数（2)的形式十分复杂，本实施例中仅考虑Si =S#情况。因此，根据变量kJPk2的值可能有下述情形：
[0080] 情形 1 %= 0,k2= 0):
[0081] 在这种情况下，m=M/Si，n=R/Sj，式（2)的评估函数可简化为式（2-5)。
[0083] 式（2-5)中，\3以，SJ表示子矩阵①乘子矩阵③的总计算时钟节拍数，M表示矩 阵A的行数，N表示矩阵A的列数，R表示矩阵B的列数，5,表示矩阵A被按行划分的行数。
[0084] 情形 2〇^=0,k2=l):
[0085] 在这种情况下，m=M/Sp 式⑵的评估函数可简化为式（2- 6)。
[0086]
[0087] 式（2-6)中，表示子矩阵①乘子矩阵③的总计算时钟节拍数， USpSp表示子矩阵①乘子矩阵④的总计算时钟节拍数，M表示矩阵A的行数，N表示矩 阵A的列数，R表示矩阵B的列数，S」表示矩阵B被按列划分的列数。
[0088] 情形 3 %= 1，k 2= 0):
[0089] 在这种情况下，与情形2类似，m= \jn=R/Sj，式（2)的评估函数可简化 为式（2-7)。
[0091] 式（2-7)中，TjSpSj表示子矩阵①乘子矩阵③的总计算时钟节拍数， T2,3(Si，Sp表示子矩阵②乘子矩阵③的总计算时钟节拍数，M表示矩阵A的行数，N表示矩 阵A的列数，R表示矩阵B的列数，5,表示矩阵A被按行划分的行数。
[0092] 情形 4(k1= l，k2= 1):
[0093] 令S'i=M-SiXm，S'」=R-SjXn，为了更好地简化讨论并且不失一般性，假设 S'i>S'j，式（2)的评估函数可简化为式（2-8)。
[0095] 式（2-8)中，USpSj)表示子矩阵①乘子矩阵③的总计算时钟节拍数， USi，sp表示子矩阵①乘子矩阵④的总计算时钟节拍数，T2,3(Si，sp表示子矩阵②乘子矩 阵③的总计算时钟节拍数，T2,4(Si，SP表示子矩阵②乘子矩阵④的总计算时钟节拍数。
[0096] 对于情形1至3,当Si的值小于或等于处理单元链中处理单元的个数时，评估函数 可取得最小值，因此在这三种情况下，最优的分块大小Si应当小于或等于处理单元链中处 理单元的个数P。对于情形4,无法直接求得评估函数的最小值，因此必须借助数学工具对 其进行分析，求得最优分块。对于具体的矩阵规模，M，N，R的值是已知的，也就意味着评估 函数是一元函数。本实施例中，具体是使用MATLAB数学工具对评估函数生成满足式（1)约 束的可视化图像，以矩阵A被按行划分的行数Si和矩阵B被按列划分的列数S/目等为前 提，可以发现评估函数的自变量（分块大小）只能取整数值，并且在整数点函数导数为零， 说明评估函数存在很多局部最值，本实施例通过作图分析的方式确定最优的分块大小 通过可视化图像确定评估函数定义域内的全局最小值，得到分块大小集合内最优的分块大 小Si〇
[0097] 以规模为M= 128、N= 1728、R= 169的矩阵乘为例，此时评估函数的表达式具体 如式（2-9)所示，并且自变量Si满足约束条件如式（2-10)所示：
[0099]
[0100] 式（2-9)中，51表示矩阵A被按行划分的行数，即最优的分块大小；L128/&」表示 对128/Si的结果向下取整；max为取较大值函数。本实施例中，具体是使用MATLAB数学工 具对评估函数生成满足式（1)约束的可视化图像如图4所示。从图4可知，最优的分块大 小Si在取值43时，评估函数达到最小值。子块大小Si的选取对矩阵加速器的性能发挥影 响很大，只有唯一的一个最优的分块大小使矩阵加速器发挥最高计算效率。
[0101] 本实施例中，步骤3)的详细步骤包括：
[0102] 3. 1)将矩阵A以规模为SjN的子块为单位按行分块，得到包含整数个规模为SjN 的子块的子矩阵①，如果仍有剩余的不规则的子块，则将剩余的不规则的子块作为子矩阵 ②；将矩阵B以规模为N*Sd^子块为单位按行分块，得到包含整数个规模为N*S,的子块的 子矩阵③，如果仍有剩余的不规则的子块，则将剩余的不规则的子块作为子矩阵④；将矩阵 A和矩阵B的矩阵乘运算视为得到的所有子矩阵的两两相乘，根据得到的子矩阵数量建立 使用指针相连接的DMA描述符链表，使得每一对相乘的子矩阵对应一个DMA描述符链表；
[0103] 3. 2)选择一个子块乘运算作为当前子块乘运算；
[0104] 3. 3)获取当前子块乘运算在矩阵A中对应子块数据的首地址、传输长度及传输步 长，获取当前子块乘运算在矩阵B中对应子块数据的首地址、传输长度及传输步长，将当前 子块乘运算在矩阵A中对应子块数据的首地址、传输长度及传输步长和当前子块乘运算在 矩阵B中对应子块数据的首地址、传输长度及传输步长封装生成一个DMA描述符，根据当前 子块乘运算所属的一对相乘的子矩阵确定对应的DMA描述符链表，并将该DMA描述符写入 对应的DMA描述符链表中；
[0105] 3. 4)判断是否所有子块乘运算已经完成处理，如果尚未完成所有子块乘运算的处 理，则选择下一个子块乘运算作为当前子块乘运算，跳转执行步骤3. 3);否则，如果已经完 成所有子块乘运算的处理则跳转执行步骤4)。
[0106] 假设矩阵A和矩阵B是两个规模为4*4的矩阵，当最优的子块大小Si= 2时，将矩 阵A分块为包含整数个规模为S#N的子块的子矩阵①、没有剩余的不规则的子矩阵②，将 矩阵B分块为包含整数个规模为N*Si的子块的子矩阵③、没有剩余的不规则的子矩阵④， 因此只存在一对相乘的子矩阵子矩阵①乘子矩阵③，此时矩阵共被划分为四个子块，矩阵 乘A*B总共需要计算四次分块乘法（A与B子块两两相乘）。每对子块乘法对应一个DMA描 述符（BufferDescriptor,BD)，因此一共有BD0~BD4共四个BD，四个BD组织成一条DMA 描述符链表存于主存中，如图5所示。DMA根据当前BD的信息就能够知道下一个BD在内 存的位置，因此DMA只需知道第一个BD的首地址即可完成对所有BD的读取，每个BD包含 了对应子块乘法的操作数的读取信息，DMA正是根据这些信息读取相应的数据。在整个读 数过程中，CPU无需参与，只需等待最终结果写回即可，这极大减少了主机和DMA的通信开 销。需要说明的是，如果出现子块数不是整数时，我们需要配置多个BD链，最差的情况下需 要配置四条BD链：最差的情况下，将矩阵A分块为包含整数个规模为SjN的子块的子矩阵 ①和剩余的不规则的子矩阵②、将矩阵B分块为包含整数个规模为N*Si的子块的子矩阵③ 和剩余的不规则的子矩阵④，因此共包含子矩阵①乘子矩阵③、子矩阵①乘子矩阵④、子矩 阵②乘子矩阵③、子矩阵②乘子矩阵④共四对相乘的子矩阵，此时则需要生成四个DMA描 述符链表，相邻DMA描述符链表之间使用指针相连接的DMA描述符链表。
[0107] 本实施例中，步骤4)的详细步骤包括：
[0108] 4. 1)CPU将DMA描述符链表中第一个DMA描述符的首地址配置给矩阵乘加速器的 DMA，矩阵乘加速器的DMA开始根据首地址读取第一个DMA描述符作为 当前描述符；
[0109] 4. 2)矩阵乘加速器的DMA解析当前描述符，得到当前描述符对应子块乘运算所需 的一对子块在外存中的地址及传输长度，并根据预设的传输步长和当前描述符中携带的一 对子块的地址及传输长度，以分时的方式交换读取当前描述符对应子块乘运算所需数据， 并将输入数据存入FIFO缓存；
[0110] 4.3)基于FIFO缓存中当前描述符对应子块乘运算所需数据，通过矩阵乘加速器 中至少一条处理单元链中的前Si个处理单元进行当前描述符对应子块乘运算；
[0111] 4. 4)将当前描述符对应子块乘运算的结果通过矩阵乘加速器的DMA分别写回主 存；
[0112] 4. 5)判断DMA描述符链表中的所有DMA描述符是否已经处理完毕，如果尚未处理 完毕，则从DMA描述符链表中选择下一个DMA描述符作为当前描述符，跳转执行步骤4. 2); 如果已经处理完毕，则结束并退出。
[0113] 本实施例中，步骤4. 3)的详细步骤包括：
[0114] 4. 3. 1)在矩阵乘加速器中确定一条用于当前描述符对应子块乘运算的处理单元 链，处理单元链至少包含Si个处理单元；矩阵乘加速器中设有用于控制各个处理单元工作 状态的状态机，状态机包括预取阶段、预取及计算阶段、提交阶段共三个状态；初始化变量 k为0,控制状态机进入预取状态；
[0115] 4. 3. 2)在预取阶段，针对当前描述符对应子块乘运算，矩阵乘加速器的DMA将属 于矩阵A的子块中的第k列数据附带上Sif处理单元的编号信息后发出给处理单元链，处 理单元链的前Si个处理单元根据编号信息预取编号对应的数据，控制状态机进入预取及计 算阶段；
[0116] 4. 3. 3)在预取及计算阶段，矩阵乘加速器的DMA将属于矩阵B的子块中第k行数 据以数据流的形式依次通过处理单元链的前Si个处理单元，且前Sif处理单元每一拍将接 收到的数据和属于矩阵A的子块中的第k列数据部分进行乘法运算，并将乘法运算结果和 片上存储中存储的上一次累加结果进行累加，将累加结果写回片上存储作为下一次累加的 输入；同时，处理单元链的前Sif处理单元分别根据编号信息预取属于矩阵A的子块中的 第k+1列数据；
[0117] 4. 3. 4)判断变量k的值是否等于N减1，如果变量k的值等于N减1，则控制状态 机进入提交阶段，跳转执行步骤4. 3. 5);否则，将初始化变量k加1，控制状态机进入预取及 计算阶段，跳转执行步骤4. 3. 3);
[0118] 4. 3. 5)在提交阶段下，处理单元链的前Si个处理单元分别将最后一次累加的结果 递给矩阵乘加速器的DMA，最终由矩阵乘加速器的DMA将结果写入主存中的指定区域。
[0119] 如图6所示，本实施例在步骤3)中，CPU根据矩阵信息生成各个子块乘运算对应 的DMA描述符（BufferDescriptor,BD)，BD组织成四条DMA描述符链表存于主存中；DMA 根据当前BD的信息就能够知道下一个BD在内存的位置，因此DMA只需知道第一个BD的首 地址即可完成对所有BD的读取，每个BD包含了对应子块乘法的操作数的读取信息，DMA正 是根据这些信息读取相应的数据。在整个读数过程中，CPU无需参与，只需等待最终结果写 回即可，这极大减少了主机和DMA的通信开销。在进入步骤4)后，首先CPU将BDO(第一条 DMA描述符链表的第一个DMA描述符）的首地址配置给DMA并启动DMA。DMA在工作状态 下，通过数据流的形式，自动读取DMA描述符链表中的DMA描述符至FIFO中，当读到第四个 DMA描述符链表BD4的最后一个DMA描述符时，DMA读取结束，等待矩阵加速器将结果写回 主存。
[0120] 如图7所示，本实施例中包含矩阵加速器的加速器系统基于xilinx公司的Zynq XC7Z045芯片上实现，主要包括处理机系统、片外存储和矩阵乘加速器，其中处理机系统包 括CPU、AXI总线模块和存储控制器，CPU具体采用ARM双核A9处理器，CPU可以通过AXI 总线模块经由存储控制器访问主存。处理机系统为矩阵加速器提供了AXI总线高速访存接 口，矩阵加速器通过该接口也能访问主存，并与CPU共享主存空间。矩阵加速器由DMA和处 理单元（ProcessingElement，PE)链（PE0~PEn)组成，每一个DMA和一个处理单元链相 连，且由链首的PE0直接和DMA交互，DMA直接与高速访存接口相连接，可直接从主存读取 加速器所需的数据。处理单元链是由若干个处理单元（PE)组成的链式结构，每个处理单元 的结构完全一致，数据在相邻的处理单元间进行传递。初始化的时候处理数据存储在片外 存储中，计算过程中的中间数据存在片上存储中，最终的结果又写回片外存储。具体过程是 第一个处理单元PE0从DMA读取数据，依次向后一个PE传递数据，每个PE保存中间结果， 等到全部计算完毕后，计算结果以相反方向向前一个PE传递，并最终由DMA将结果写回主 存。矩阵加速器可扩展为若干条处理单元链，每条处理单元链都有其对应的DMA。由于高性 能接口之间相对独立，都能够提供稳定的带宽，因此在DMA协助下，处理单元链之间可以高 效并行地完成计算任务。
[0121] DMA主要用以管理片外存储于加速器间的数据传输。它的主要特点是计算大规模 矩阵时往往只需要CPU配置一次，配置信息足够完整并存于片外存储中，DMA自动读取配置 信息后就能够完成所有的数据传输工作。参见图6,DMA的工作步骤主要包括：步骤1，CPU 根据矩阵分块后的矩阵的数据首地址、传输长度等信息生成DMA描述符，并将其组织成为 链表结构，存入片外存储。一个矩阵子块乘法对应一个描述符，描述符之间用指针相连接。 步骤2,CPU将描述符链表的在外存的首地址配置给DMA。DMA根据该地址自动读取DMA描 述符，然后解析出描述符的配置信息，得到矩阵输入数据在外存的地址，并根据传输步长， 传输长度等信息，以分时方式交换读取矩阵的输入数据。获得输入数据后，DMA将其存入相 应的FIFO缓存，为加速器提供操作数据。步骤3,DMA读取完一个子块乘法所需的操作数之 后，可以根据下个描述符的地址（当前描述符的信息）继续读取下一个子块乘法的操作数， 在此过程中，DMA还负责将上一个子块乘法的运算结果写回外存。步骤4,重复步骤3直到 DMA读取到最后一个DMA描述符。
[0122] 如图8所示，本实施例中每个处理单元由计算模块、计算数据传送模块、数据存储 模块和阶段同步控制模块组成。
[0123] 计算单元采用Xilinx公司提供的浮点运算IP核进行实现。计算单元包括乘法器 和加法器，其中加法器的流水线级别为11。
[0124] 计算数据传送模块由FIF0_A、FIF0_B、FIF0_C共三个FIFO队列组成，FIF0_A和 FIFO_B负责将输入数据传递到下一个PE。FIFO_C负责将结果传递到前一个PE中。
[0125] 数据存储模块主要包括一个双端口BRAM(MEM_C)和地址生成器。当乘法器产生第 一个结果时，地址生成器产生读信号和读地址，从MEM_C读取加法器的另一个操作数，也就 是上一阶段的临时数据送往加法器；当加法器的第一个结果产生时，地址生成器产生写信 号和写地址，结果被再次存入MEM_C。
[0126] 阶段同步控制模块内部主要实现了两个计数器，当矩阵A子块的新的一列和矩阵 B子块新的一行进入PE时，两个计数器分别开始计数。每个计数器都有相应的阈值，一旦 其中一个计数器达到阈值，阶段同步控制模块就会使当前PE停止从前一个PE读取相对应 的数据。这样一来达到阈值的计数器就会停止计数，一直到另一个计数器也达到其阈值， 两个计数器同时被清零，重新开始为下一阶段的数据计数。通过阶段同步保证了矩阵A的 某一列的第一个数据总是和矩阵B对应行的第一个数据同时进入每个PE。每个PE都有一 个ID(PID)，矩阵A的每个数据在进入矩阵链之前被加入了编号信息用以表征这个数据归 属于哪个处理单元。阶段同步控制模块，用以解决支持可变分块引起的数据不同步问题；此 外，每个处理单元还包括用以控制与相邻处理单元间的数据交互的控制逻辑。
[0127] 需要说明的是，矩阵加速器中处理单元的结构实现并不局限于如图8所示的特定 结构，毫无疑问，本领域技术人员也可以根据需要采用其他结构的处理单元来实现矩阵加 速器中的子块乘法运算，故在此不再赘述。
[0128] 根据图2的原理可知，编号信息依次为0，1，…，Sri.iS'P时，那些满足PID> =PE不允许启动工作。我们通过对数据编号和PID进行对比，使得当数据编号小于 当前PID时，数据才允许被写入FIFO_A并传递到下一个PE。这样一来，数据就只传递到前 SifPE中，我们的矩阵乘法器就能支持可变大小的分块。
[0129] 本实施例中，每个处理单元内部实现了 一个状态机用来控制数据流。处理单元工 作时主要有三个阶段，我们以处理单元链中的第一个处理单元（PE0)为例：第一，预取阶 段。PE0从DMA的FIFO中读取矩阵A的第一列数据，并通过FIFO_A传递给相邻的PE。数 据从FIF0_A出来后在高位加入了编号信息，每个PE根据自己的ID号预取各自的数据。第 二，计算-预取阶段。PE0开始同时读取矩阵A的第k+1列以及矩阵B的第k行数据并通过 FIF0_A和FIF0_B传递给相邻的PE;矩阵B的数据与预取的A数据相乘，同时更新A的另一 个预取缓冲（双缓冲，TA0，TA1)。经过一定延迟后对乘法结果进行加法运算，其中另一个操 作数来自MEM_C，加法计算完毕后，又将结果存入MEM_C。第三阶段，提交阶段。每个PE将 最后一次累加结果存入FIF0_C，同时MEM_C内容清零。存入完毕后 ，每个PE从后一个PE读 取结果存入FIF0_C，由此结果数据依次前移至DMA的FIF0_C中。
[0130] 本实施例以典型卷积神经网络（CNN)中如表1所示的五个非均匀矩阵规模 Convl~Conv5为例，应用本实施例以及经典分块算法分别对表1中的五个规模的矩阵乘进 行加速对比，表2是本实施例中每个卷积层的矩阵规模通过优化分块策略算出的最优的分 块大小Sp
[0131] 表1 :五个非均匀矩阵规模实例表。
[0133] 表2 :五个非均匀矩阵规模实例的最优的分块大小31数据表。
[0135] 最终，得到的性能对比图如图9所示。参见图9可知，当矩阵规模不均匀程度比较 高时（矩阵的行列大小相差较大），本实施例的矩阵乘甲酸方法能达到很高的性能提升，例 如对于Conv5,本实施例的矩阵乘甲酸方法相对于经典分块算法有大约12%的性能提升。
[0136] 以上所述仅是本发明的优选实施方式，本发明的保护范围并不仅局限于上述实施 例，凡属于本发明思路下的技术方案均属于本发明的保护范围。应当指出，对于本技术领域 的普通技术人员来说，在不脱离本发明原理前提下的若干改进和润饰，这些改进和润饰也 应视为本发明的保护范围。
【主权项】
1. 一种支持可变分块的矩阵乘加速方法，其特征在于步骤包括： 1) 输入矩阵乘运算所需的M*N的矩阵A和N*R的矩阵B ; 2) 根据矩阵A和矩阵B的规模确定子块大小Si,将矩阵A以规模为SjN的子块为单位 进行按行分块，将矩阵B以规模为N*Si的子块为单位进行按列分块，使得矩阵乘运算等同 为多个子块乘运算； 3) 为每一个子块乘运算所需数据生成一个DM描述符，将所有子块乘运算的DM描述 符构建DM描述符链表并存入主存； 4) 针对每一个子块乘运算，通过矩阵乘加速器的DM从主存的DM描述符链表读取子 块乘运算所需数据，然后通过矩阵乘加速器中至少一条处理单元链中的前Si个处理单元进 行子块乘运算，并将各个子块乘运算的结果通过矩阵乘加速器的DM分别写回主存。2. 根据权利要求1所述的支持可变分块的矩阵乘加速方法，其特征在于，所述步骤2) 中根据矩阵A和矩阵B的规模确定子块大小详细步骤包括： 2. 1)输入矩阵A和矩阵B的矩阵规模，所述矩阵规模包括矩阵A的行数M、矩阵A的列 数N、矩阵B的列数R三者的值； 2. 2)根据矩阵规模获取满足式（1)所示约束的分块大小集合；式⑴中，示矩阵A被按行划分的行数，S廣示矩阵B被按列划分的列数，M表 示矩阵A的行数，N表示矩阵A的列数，BW表示矩阵乘加速器中单条处理单元链的访 存带宽，HiaxiSi, SjI表示从SjP S」中取较大值，F表示访存频率；P表示矩阵乘加速器 中单条处理单元链的处理单元数量，Stageadd表示矩阵乘加速器中的加法器流水级数， Hiaxm-SiXm, R-SjXn}表示从M-SiXm和R-SjXn中取较大值，R表示矩阵B的列数，m表示 矩阵A的行数M除以矩阵A被按行划分的行数结果的向下取整值，η表示矩阵B的列 数R除以矩阵B被按列划分的列数\_的结果的向下取整值； 2. 3)将矩阵A视为包含整数个规模为S#N的子块的子矩阵①和剩余的不规则的子矩 阵②、将矩阵B视为包含整数个规模为N*Sj^子块的子矩阵③和剩余的不规则的子矩阵 ④，将矩阵A和矩阵B的矩阵乘运算视为子矩阵①~子矩阵④四者中的两两相乘，建立式 (2)所示的评估函数； f (Si, Sj) = 8,.)+^3(8^ Sj) (2) 式（2)中，f (Si, Sp表示评估函数的值，51表示矩阵A被按行划分的行数，S ^表示矩阵 B被按列划分的列数，表示子矩阵①乘子矩阵③的计算时钟节拍数，Tli4(SilSj) 表示子矩阵①乘子矩阵④的计算时钟节拍数，T2i3(Si，sp表示子矩阵②乘子矩阵③的计算 时钟节拍数，(Si, sp表示子矩阵②乘子矩阵④的计算时钟节拍数； 2.4)使用MATLAB数学工具对式（2)所示的评估函数生成满足式（1)约束的可视化图 像，以矩阵A被按行划分的行数Si和矩阵B被按列划分的列数S d目等为前提，通过所述可 视化图像确定评估函数定义域内的全局最小值，得到分块大小集合内最优的分块大小Sp3.根据权利要求2所述的支持可变分块的矩阵乘加速方法，其特征在于，所述步骤3) 的详细步骤包括： 3. 1)将矩阵A以规模为S#N的子块为单位按行分块，得到包含整数个规模为S#N的 子块的子矩阵①，如果仍有剩余的不规则的子块，则将剩余的不规则的子块作为子矩阵②； 将矩阵B以规模为N*Sj^子块为单位按行分块，得到包含整数个规模为N*S i的子块的子矩 阵③，如果仍有剩余的不规则的子块，则将剩余的不规则的子块作为子矩阵④；将矩阵A和 矩阵B的矩阵乘运算视为得到的所有子矩阵的两两相乘，根据得到的子矩阵数量建立使用 指针相连接的DM描述符链表，使得每一对相乘的子矩阵对应一个DM描述符链表； 3. 2)选择一个子块乘运算作为当前子块乘运算； 3. 3)获取当前子块乘运算在矩阵A中对应子块数据的首地址、传输长度及传输步长， 获取当前子块乘运算在矩阵B中对应子块数据的首地址、传输长度及传输步长，将当前子 块乘运算在矩阵A中对应子块数据的首地址、传输长度及传输步长和当前子块乘运算在矩 阵B中对应子块数据的首地址、传输长度及传输步长封装生成一个DM描述符，根据当前子 块乘运算所属的一对相乘的子矩阵确定对应的DM描述符链表，并将该DM描述符写入对 应的DM描述符链表中； 3. 4)判断是否所有子块乘运算已经完成处理，如果尚未完成所有子块乘运算的处理， 则选择下一个子块乘运算作为当前子块乘运算，跳转执行步骤3. 3);否则，如果已经完成 所有子块乘运算的处理则跳转执行步骤4)。4. 根据权利要求3所述的支持可变分块的矩阵乘加速方法，其特征在于，所述步骤4) 的详细步骤包括： 4. I) CPU将DM描述符链表中第一个DM描述符的首地址配置给矩阵乘加速器的DMA， 矩阵乘加速器的DM开始根据首地址读取第一个DM描述符作为当前描述符； 4. 2)矩阵乘加速器的DM解析当前描述符，得到当前描述符对应子块乘运算所需的一 对子块在外存中的地址及传输长度，并根据预设的传输步长和当前描述符中携带的一对子 块的地址及传输长度，以分时的方式交换读取当前描述符对应子块乘运算所需数据，并将 输入数据存入FIFO缓存； 4. 3)基于FIFO缓存中当前描述符对应子块乘运算所需数据，通过矩阵乘加速器中至 少一条处理单元链中的前Si个处理单元进行当前描述符对应子块乘运算； 4. 4)将当前描述符对应子块乘运算的结果通过矩阵乘加速器的DM分别写回主存； 4. 5)判断DMA描述符链表中的所有DMA描述符是否已经处理完毕，如果尚未处理完毕， 则从DM描述符链表中选择下一个DM描述符作为当前描述符，跳转执行步骤4. 2);如果 已经处理完毕，则结束并退出。5. 根据权利要求4所述的支持可变分块的矩阵乘加速方法，其特征在于，所述步骤 4. 3)的详细步骤包括： 4. 3. 1)在矩阵乘加速器中确定一条用于当前描述符对应子块乘运算的处理单元链，所 述处理单元链至少包含Si个处理单元；所述矩阵乘加速器中设有用于控制各个处理单元工 作状态的状态机，所述状态机包括预取阶段、预取及计算阶段、提交阶段共三个状态；初始 化变量k为0,控制状态机进入预取状态； 4. 3. 2)在预取阶段，针对当前描述符对应子块乘运算，矩阵乘加速器的DM将属于矩 阵A的子块中的第k列数据附带上SiA处理单元的编号信息后发出给处理单元链，处理单 元链的前Si个处理单元根据编号信息预取编号对应的数据，控制状态机进入预取及计算阶 段； 4. 3. 3)在预取及计算阶段，矩阵乘加速器的DM将属于矩阵B的子块中第k行数据以 数据流的形式依次通过处理单元链的前Si个处理单元，且前S i个处理单元每一拍将接收 到的数据和属于矩阵A的子块中的第k列数据部分进行乘法运算，并将乘法运算结果和片 上存储中存储的上一次累加结果进行累加，将累加结果写回片上存储作为下一次累加的输 入；同时，处理单元链的前SiA处理单元分别根据编号信息预取属于矩阵A的子块中的第 k+Ι列数据； 4. 3. 4)判断变量k的值是否等于N减1，如果变量k的值等于N减1，则控制状态机进 入提交阶段，跳转执行步骤4. 3. 5);否则，将初始化变量k加1，控制状态机进入预取及计算 阶段，跳转执行步骤4.3.3); 4. 3. 5)在提交阶段下，处理单元链的前SiA处理单元分别将最后一次累加的结果递给 矩阵乘加速器的DMA，最终由矩阵乘加速器的DMA将结果写入主存中的指定区域。
【专利摘要】本发明公开了一种支持可变分块的矩阵乘加速方法，步骤包括：输入矩阵A和矩阵B，根据矩阵A和矩阵B的规模确定子块大小Si，将矩阵A以规模为Si*N的子块为单位进行按行分块，将矩阵B以规模为N*Si的子块为单位进行按列分块，为每一个子块乘运算所需数据生成一个DMA描述符，将所有DMA描述符构建DMA描述符链表；针对每一个子块乘运算，根据主存的DMA描述符链表读取子块乘运算所需数据，通过矩阵乘加速器中的一条处理单元链进行子块乘运算，并将结果通过DMA写回主存。本发明具有支持可变分块，可根据分块大小调整使用的处理单元数目，加速非均匀矩阵乘运算加速效率高的优点。
【IPC分类】G06F7/523, G06F17/16
【公开号】CN104899182
【申请号】CN201510312188
【发明人】文梅, 沈俊忠, 乔寓然, 杨乾明, 苏华友, 肖涛, 陈照云, 张春元
【申请人】中国人民解放军国防科学技术大学
【公开日】2015年9月9日
【申请日】2015年6月9日
最新回复(0)