CN104899182B - 一种支持可变分块的矩阵乘加速方法 - Google Patents

Abstract

本发明公开了一种支持可变分块的矩阵乘加速方法，步骤包括：输入矩阵A和矩阵B，根据矩阵A和矩阵B的规模确定矩阵A被按行划分的行数S_i，将矩阵A以规模为S_i*N的子块为单位进行按行分块，将矩阵B以规模为N*S_i的子块为单位进行按列分块，为每一个子块乘运算所需数据生成一个DMA描述符，将所有DMA描述符构建DMA描述符链表；针对每一个子块乘运算，根据主存的DMA描述符链表读取子块乘运算所需数据，通过矩阵乘加速器中的一条处理单元链进行子块乘运算，并将结果通过DMA写回主存。本发明具有支持可变分块，可根据分块大小调整使用的处理单元数目，加速非均匀矩阵乘运算加速效率高的优点。

Description

一种支持可变分块的矩阵乘加速方法

技术领域

本发明涉及嵌入式平台下的矩阵乘加速技术，具体涉及一种支持可变分块的矩阵乘加速方法。

背景技术

随着半导体制造工艺的发展以及集成电路技术的进步，单芯片上能够集成越来越多的晶体管，使用可编程器件，特别是FPGA(Field Programmable Gate Array)芯片进行设计成为了目前构建嵌入式系统以及硬件加速平台的一种重要方式。当前FPGA芯片提供了专用的算术模块、大量的逻辑资源和存储资源，以及外部存储器接口、网络接口和其它外围接口，为构建高性能计算系统提供了条件，也使FPGA可重构计算系统成为加速科学计算的一种重要选择。当前的FPGA芯片相对于DSP(Digital Signal Processor,数字信号处理器)具有可编程的优势，同时能够并行处理海量数据，既具有通用处理器的灵活性，又具有ASIC(Application Specific Integrated Circuit,专用集成电路)的高性能，在嵌入式计算领域备受青睐。

浮点矩阵乘法是数字信号的基本算法，同时也是许多科学计算方法的基本运算。在数字图像处理，计算机视觉的快速处理以及工业实时控制等领域都被广泛应用。但由于实际应用中，浮点矩阵规模通常较大，矩阵乘算法本身复杂度较高、处理效率较低，成为限制系统性能提升的瓶颈所在，因此为此类应用设计高性能的硬件结构是当前FPGA结构设计的研究热点。

尽管近年来提出了很多基于FPGA的矩阵乘加速器设计，但是都缺乏对非均匀的大规模矩阵加速的讨论和支持，这种大规模矩阵的特征是行列数相差很大(>＝10倍)，并且广泛存在于很多现代应用领域中，如图像处理，深度学习等。在这些应用领域之中，矩阵乘占据了计算量的主要部分。由于单片FPGA芯片上存储资源和计算资源十分有限，加速大规模矩阵乘时往往需要对矩阵进行分块。对于链式结构的矩阵乘加速器，虽然对大多数大规模矩阵的加速效果十分明显，但对于加速非均匀矩阵时的计算效率却很低，其主要原因就是这类加速器往往只支持固定的分块，也就是说，分块大小与矩阵链长(矩阵链中处理单元个数)相等或者是其倍数。当分块大小与矩阵加速器链长不匹配时，加速器计算效率会发生明显的下降。据我们所知，至今还没有公开文献涉及支持可变分块的矩阵乘法器设计，也没有经典分块算法基础上的分块优化问题方面的相关研究，因此如何选择最优的分块，使得矩阵乘加速器获得更高的计算效率，以更好地适应现代应用的加速需求，已经成为一项亟待解决的关键技术问题。

发明内容

本发明要解决的技术问题是：针对现有技术的上述技术问题，提供一种支持可变分块，可根据分块大小调整使用的处理单元数目，加速非均匀矩阵乘运算加速效率高的支持可变分块的矩阵乘加速方法。

为了解决上述技术问题，本发明采用的技术方案为：

一种支持可变分块的矩阵乘加速方法，步骤包括：

1)输入矩阵乘运算所需的M*N的矩阵A和N*R的矩阵B；

2)根据矩阵A和矩阵B的规模确定矩阵A被按行划分的行数S_i，将矩阵A以规模为S_i*N的子块为单位进行按行分块，将矩阵B以规模为N*S_i的子块为单位进行按列分块，使得矩阵乘运算等同为多个子块乘运算；

3)为每一个子块乘运算所需数据生成一个DMA描述符，将所有子块乘运算的DMA描述符构建DMA描述符链表并存入主存；

4)针对每一个子块乘运算，通过矩阵乘加速器的DMA从主存的DMA描述符链表读取子块乘运算所需数据，然后通过矩阵乘加速器中至少一条处理单元链中的前S_i个处理单元进行子块乘运算，并将各个子块乘运算的结果通过矩阵乘加速器的DMA分别写回主存。

优选地，所述步骤2)中根据矩阵A和矩阵B的规模确定矩阵A被按行划分的行数S_i的详细步骤包括：

2.1)输入矩阵A和矩阵B的矩阵规模，所述矩阵规模包括矩阵A的行数M、矩阵A的列数N、矩阵B的列数R三者的值；

2.2)根据矩阵规模获取满足式(1)所示约束的分块大小集合；

式(1)中，S_i表示矩阵A被按行划分的行数，S_j表示矩阵B被按列划分的列数，M表示矩阵A的行数，N表示矩阵A的列数，BW表示矩阵乘加速器中单条处理单元链的访存带宽，max{S_i,S_j}表示从S_i和S_j中取较大值，F表示访存频率；P表示矩阵乘加速器中单条处理单元链的处理单元数量，Stage_add表示矩阵乘加速器中的加法器流水级数，max{M-S_i×m,R-S_j×n}表示从M-S_i×m和R-S_j×n中取较大值，R表示矩阵B的列数，m表示矩阵A的行数M除以矩阵A被按行划分的行数S_i的结果的向下取整值，n表示矩阵B的列数R除以矩阵B被按列划分的列数S_j的结果的向下取整值；

2.3)将矩阵A视为包含整数个规模为S_i*N的子块的子矩阵①和剩余的不规则的子矩阵②、将矩阵B视为包含整数个规模为N*S_i的子块的子矩阵③和剩余的不规则的子矩阵④，将矩阵A和矩阵B的矩阵乘运算视为子矩阵①～子矩阵④四者中的两两相乘，建立式(2)所示的评估函数；

f(S_i,S_j)＝T_1,3(S_i,S_j)+T_1,4(S_i,S_j)+T_2,3(S_i,S_j)+T_2,4(S_i,S_j) (2)

式(2)中，f(S_i,S_j)表示评估函数的值，S_i表示矩阵A被按行划分的行数，S_j表示矩阵B被按列划分的列数，T_1,3(S_i,S_j)表示子矩阵①乘子矩阵③的计算时钟节拍数，T_1,4(S_i,S_j)表示子矩阵①乘子矩阵④的计算时钟节拍数，T_2,3(S_i,S_j)表示子矩阵②乘子矩阵③的计算时钟节拍数，T_2,4(S_i,S_j)表示子矩阵②乘子矩阵④的计算时钟节拍数；

2.4)使用MATLAB数学工具对式(2)所示的评估函数生成满足式(1)约束的可视化图像，以矩阵A被按行划分的行数S_i和矩阵B被按列划分的列数S_j相等为前提，通过所述可视化图像确定评估函数定义域内的全局最小值，得到分块大小集合内最优的矩阵A被按行划分的行数S_i。

优选地，所述步骤3)的详细步骤包括：

3.1)将矩阵A以规模为S_i*N的子块为单位按行分块，得到包含整数个规模为S_i*N的子块的子矩阵①，如果仍有剩余的不规则的子块，则将剩余的不规则的子块作为子矩阵②；将矩阵B以规模为N*S_i的子块为单位按行分块，得到包含整数个规模为N*S_i的子块的子矩阵③，如果仍有剩余的不规则的子块，则将剩余的不规则的子块作为子矩阵④；将矩阵A和矩阵B的矩阵乘运算视为得到的所有子矩阵的两两相乘，根据得到的子矩阵数量建立使用指针相连接的DMA描述符链表，使得每一对相乘的子矩阵对应一个DMA描述符链表；

3.2)选择一个子块乘运算作为当前子块乘运算；

3.3)获取当前子块乘运算在矩阵A中对应子块数据的首地址、传输长度及传输步长，获取当前子块乘运算在矩阵B中对应子块数据的首地址、传输长度及传输步长，将当前子块乘运算在矩阵A中对应子块数据的首地址、传输长度及传输步长和当前子块乘运算在矩阵B中对应子块数据的首地址、传输长度及传输步长封装生成一个DMA描述符，根据当前子块乘运算所属的一对相乘的子矩阵确定对应的DMA描述符链表，并将该DMA描述符写入对应的DMA描述符链表中；

3.4)判断是否所有子块乘运算已经完成处理，如果尚未完成所有子块乘运算的处理，则选择下一个子块乘运算作为当前子块乘运算，跳转执行步骤3.3)；否则，如果已经完成所有子块乘运算的处理则跳转执行步骤4)。

优选地，所述步骤4)的详细步骤包括：

4.1)CPU将DMA描述符链表中第一个DMA描述符的首地址配置给矩阵乘加速器的DMA，矩阵乘加速器的DMA开始根据首地址读取第一个DMA描述符作为当前描述符；

4.2)矩阵乘加速器的DMA解析当前描述符，得到当前描述符对应子块乘运算所需的一对子块在外存中的地址及传输长度，并根据预设的传输步长和当前描述符中携带的一对子块的地址及传输长度，以分时的方式交换读取当前描述符对应子块乘运算所需数据，并将输入数据存入FIFO缓存；

4.3)基于FIFO缓存中当前描述符对应子块乘运算所需数据，通过矩阵乘加速器中至少一条处理单元链中的前S_i个处理单元进行当前描述符对应子块乘运算；

4.4)将当前描述符对应子块乘运算的结果通过矩阵乘加速器的DMA分别写回主存；

4.5)判断DMA描述符链表中的所有DMA描述符是否已经处理完毕，如果尚未处理完毕，则从DMA描述符链表中选择下一个DMA描述符作为当前描述符，跳转执行步骤4.2)；如果已经处理完毕，则结束并退出。

优选地，所述步骤4.3)的详细步骤包括：

4.3.1)在矩阵乘加速器中确定一条用于当前描述符对应子块乘运算的处理单元链，所述处理单元链至少包含S_i个处理单元；所述矩阵乘加速器中设有用于控制各个处理单元工作状态的状态机，所述状态机包括预取阶段、预取及计算阶段、提交阶段共三个状态；将变量k初始化为0，控制状态机进入预取状态；

4.3.2)在预取阶段，针对当前描述符对应子块乘运算，矩阵乘加速器的DMA将属于矩阵A的子块中的第k列数据附带上S_i个处理单元的编号信息后发出给处理单元链，处理单元链的前S_i个处理单元根据编号信息预取编号对应的数据，控制状态机进入预取及计算阶段；

4.3.3)在预取及计算阶段，矩阵乘加速器的DMA将属于矩阵B的子块中第k行数据以数据流的形式依次通过处理单元链的前S_i个处理单元，且前S_i个处理单元每一拍将接收到的数据和属于矩阵A的子块中的第k列数据部分进行乘法运算，并将乘法运算结果和片上存储中存储的上一次累加结果进行累加，将累加结果写回片上存储作为下一次累加的输入；同时，处理单元链的前S_i个处理单元分别根据编号信息预取属于矩阵A的子块中的第k+1列数据；

4.3.4)判断变量k的值是否等于N减1，如果变量k的值等于N减1，则控制状态机进入提交阶段，跳转执行步骤4.3.5)；否则，将变量k加1，控制状态机进入预取及计算阶段，跳转执行步骤4.3.3)；

4.3.5)在提交阶段下，处理单元链的前S_i个处理单元分别将最后一次累加的结果递给矩阵乘加速器的DMA，最终由矩阵乘加速器的DMA将结果写入主存中的指定区域。

本发明支持可变分块的矩阵乘加速方法具有下述优点：本发明通过确定子块大小S_i，将矩阵A以规模为S_i*N的子块为单位进行按行分块，将矩阵B以规模为N*S_i的子块为单位进行按列分块，使得矩阵乘运算等同为多个子块乘运算，为每一个子块乘运算所需数据生成一个DMA描述符，将所有子块乘运算的DMA描述符构建DMA描述符链表并存入主存，针对每一个子块乘运算，通过矩阵乘加速器的DMA从主存的DMA描述符链表读取子块乘运算所需数据，然后通过矩阵乘加速器中至少一条处理单元链中的前S_i个处理单元进行子块乘运算，并将各个子块乘运算的结果通过矩阵乘加速器的DMA分别写回主存，增加了对可变分块的支持，能够根据具体的分块大小调整使用的处理单元数目，通过确定子块大小S_i来基于经典分块算法上的优化分块策略，当加速非均匀矩阵时，通过计算并应用最优分块能够使得矩阵加速器具有很高的计算效率，具有支持可变分块，可根据分块大小调整使用的处理单元数目，加速非均匀矩阵乘运算加速效率高的优点。

附图说明

图1为本发明实施例的基本流程示意图。

图2为本发明实施例将进行矩阵乘运算等同为多个子块乘运算的原理示意图。

图3为本发明实施例将进行矩阵乘运算划分为子矩阵相乘的原理示意图。

图4为本发明实施例中使用MATLAB数学工具生成的评估函数可视化图像。

图5为本发明实施例中生成的一个DMA描述符链表的结构示意图。

图6为本发明实施例中矩阵加速器的DMA的工作流程示意图。

图7为本发明实施例中加速器系统的框架结构示意图。

图8为本发明实施例中处理单元的数据流向结构示意图。

图9为本发明实施例和采用经典分块算法的性能对比图。

具体实施方式

如图1所示，本实施例支持可变分块的矩阵乘加速方法的步骤包括：

1)输入矩阵乘运算所需的M*N的矩阵A和N*R的矩阵B；

2)根据矩阵A和矩阵B的规模确定矩阵A被按行划分的行数S_i，将矩阵A以规模为S_i*N的子块为单位进行按行分块，将矩阵B以规模为N*S_i的子块为单位进行按列分块，使得矩阵乘运算等同为多个子块乘运算；

3)为每一个子块乘运算所需数据生成一个DMA描述符，将所有子块乘运算的DMA描述符构建DMA描述符链表并存入主存；

4)针对每一个子块乘运算，通过矩阵乘加速器的DMA从主存的DMA描述符链表读取子块乘运算所需数据，然后通过矩阵乘加速器中至少一条处理单元链中的前S_i个处理单元进行子块乘运算，并将各个子块乘运算的结果通过矩阵乘加速器的DMA分别写回主存。

如图2所示，对于M*N的矩阵A和N*R的矩阵B的矩阵乘运算而言，其矩阵乘运算得到的矩阵C的规模为M*R。因此，将矩阵A和N*R的矩阵B的矩阵乘运算划分为多个小矩阵相乘，其中矩阵A按行划分为多个S_i*N子块，矩阵B按列划分为多个N*S_j子块，每对小子块相乘就得到一个S_i*S_j的矩阵，也就是最后的结果C矩阵的一个子块。参见图2，现考虑子块相乘的细节，可把A的S_i*N子块看成是由N个维度为S_i的列向量组成的矩阵，B的N*S_j子块是由N个维度为S_j的行向量组成的矩阵，记第k个列向量和行向量分别为V_k和U_k(k<＝N)，把U_k和V_k相乘的结果进行累加(N次)就得到矩阵C的S_i*S_j子块。假设矩阵乘加速器有S_i个处理单元，每个处理单元每次缓冲V_k中的一个操作数，U_k中的所有操作数依次进入每个处理单元，处理单元以流水线的形式高效地进行乘加运算，并且将加法结果存入临时存储区(下一次加法运算要从存储区读取这一次的结果进行累加)，这样一来，每个处理单元并行地计算C的某一行的结果，这就是本实施例算法的核心思想。事实上当S_i＝S_j时，矩阵乘加速器进行阶段同步所导致的空拍时最少，这种情况下矩阵乘加速器的性能应该是最优的，因此本实施例仅考虑S_i＝S_j的情况。

本实施例中，步骤2)中根据矩阵A和矩阵B的规模确定矩阵A被按行划分的行数S_i的详细步骤包括：

2.1)输入矩阵A和矩阵B的矩阵规模，矩阵规模包括矩阵A的行数M、矩阵A的列数N、矩阵B的列数R三者的值；

2.2)根据矩阵规模获取满足式(1)所示约束的分块大小集合；

式(1)中，S_i表示矩阵A被按行划分的行数，S_j表示矩阵B被按列划分的列数，M表示矩阵A的行数，N表示矩阵A的列数，BW表示矩阵乘加速器中单条处理单元链的访存带宽，max{S_i,S_j}表示从S_i和S_j中取较大值，F表示访存频率；P表示矩阵乘加速器中单条处理单元链的处理单元数量，Stage_add表示矩阵乘加速器中的加法器流水级数，max{M-S_i×m,R-S_j×n}表示从M-S_i×m和R-S_j×n中取较大值，R表示矩阵B的列数，m表示矩阵A的行数M除以矩阵A被按行划分的行数S_i的结果的向下取整值，n表示矩阵B的列数R除以矩阵B被按列划分的列数S_j的结果的向下取整值；

2.3)将矩阵A视为包含整数个规模为S_i*N的子块的子矩阵①和剩余的不规则的子矩阵②、将矩阵B视为包含整数个规模为N*S_i的子块的子矩阵③和剩余的不规则的子矩阵④，将矩阵A和矩阵B的矩阵乘运算视为子矩阵①～子矩阵④四者中的两两相乘，建立式(2)所示的评估函数；

f(S_i,S_j)＝T_1,3(S_i,S_j)+T_1,4(S_i,S_j)+T_2,3(S_i,S_j)+T_2,4(S_i,S_j) (2)

式(2)中，f(S_i,S_j)表示评估函数的值，S_i表示矩阵A被按行划分的行数，S_j表示矩阵B被按列划分的列数，T_1,3(S_i,S_j)表示子矩阵①乘子矩阵③的计算时钟节拍数，T_1,4(S_i,S_j)表示子矩阵①乘子矩阵④的计算时钟节拍数，T_2,3(S_i,S_j)表示子矩阵②乘子矩阵③的计算时钟节拍数，T_2,4(S_i,S_j)表示子矩阵②乘子矩阵④的计算时钟节拍数；

2.4)使用MATLAB数学工具对式(2)所示的评估函数生成满足式(1)约束的可视化图像，以矩阵A被按行划分的行数S_i和矩阵B被按列划分的列数S_j相等为前提，通过可视化图像确定评估函数定义域内的全局最小值，得到分块大小集合内最优的矩阵A被按行划分的行数S_i。

如图3所示，将矩阵A和矩阵B的矩阵乘运算视为子矩阵①～子矩阵④四者中的两两相乘后，子矩阵①包含整数个规模为S_i*N的子块，子矩阵③包含整数个规模为N*S_i的子块，剩余的不规则的子矩阵②包含的为不规则额的子块(行数小于S_i)，剩余的不规则的子矩阵④包含的为不规则的子块(列数小于S_i)，子矩阵①～子矩阵④四者中的两两相乘即为子矩阵①乘子矩阵③、子矩阵①乘子矩阵④、子矩阵②乘子矩阵③、子矩阵②乘子矩阵④。本实施例中，由于子块乘运算要求在开始的时候预取S_i个矩阵A的数据(子块第一列)，这个过程需要S_i个节拍，每个PE计算S_i﹡S_j需要花费max{S_i,S_j}×N个节拍。由于矩阵乘加速器中的处理单元组织成链式结构，因此最后一个计算的处理单元相对于第一个处理单元有S_i个节拍的延迟，即整个处理单元流水线的建立时间。计算S_i﹡S_j子块的结果需要从主存读取(S_i+S_j)×N个数据，该过程需要时间(S_i+S_j)×N/BW，BW是单条链的访存带宽。根据所有处理单元总计算时间大于取数时间，有以下不等式(1-1)成立。

(S_i+S_j)×N/BW≤(2×S_i+max{S_i,S_j}×N)/F (1-1)

式(1-1)中，S_i表示矩阵A被按行划分的行数，S_j表示矩阵B被按列划分的列数，N表示矩阵A的列数，BW表示矩阵乘加速器中单条处理单元链的访存带宽，max{S_i,S_j}表示从S_i和S_j中取较大值，F表示访存频率。

此外，矩阵A被按行划分的行数S_i、矩阵B被按列划分的列数S_j还有以下约束：首先，根据算法，矩阵A被按行划分的行数S_i的值不能大于矩阵乘加速器中一条处理单元链中处理单元的数目；其次，考虑到流水线中的数据冲突，也就是当加法结果尚未写入片上存储器时，加法器又需要从存储器读取该结果时，加法器就会读到旧值。因此，必须同时满足式(1-2)所示的几个约束条件。

式(1-2)中，S_i表示矩阵A被按行划分的行数，P表示矩阵乘加速器中单条处理单元链的处理单元数量，max{S_i,S_j}表示从S_i和S_j中取较大值，Stage_add表示矩阵乘加速器中的加法器流水级数，max{M-S_i×m,R-S_j×n}表示从M-S_i×m和R-S_j×n中取较大值。本实施例中，Stage_add表示矩阵乘加速器中的加法器流水级数具体值为11。结合前述的式(1-1)和式(1-2)，即可推导得出本实施例前述的式(1)所示约束的分块大小集合。

假设(即m表示矩阵A的行数M除以矩阵A被按行划分的行数S_i的结果的向下取整值，n表示矩阵B的列数R除以矩阵B被按列划分的列数S_j的结果的向下取整值)，则有子矩阵①乘子矩阵③的总计算时钟节拍数T_1,3(S_i,S_j)如式(2-1)所示，子矩阵①乘子矩阵④的总计算时钟节拍数T_1,4(S_i,S_j)如式(2-2)所示，子矩阵②乘子矩阵③的总计算时钟节拍数T_2,3(S_i,S_j)如式(2-3)所示，子矩阵②乘子矩阵④的总计算时钟节拍数T_2,4(S_i,S_j)如式(2-4)所示，

T_1,3(S_i,S_j)＝m×n×(S_i+max{S_i,S_j}×N) (2-1)

T_1,4(S_i,S_j)＝k₂×m×(max{S_i,R-S_j×n}×N) (2-2)

T_2,3(S_i,S_j)＝k₁×n×(max{M-S_i×m,S_j}×N) (2-3)

T_2,4(S_i,S_j)＝k₁×k₂×(max{M-S_i×m,R-S_j×n}×N) (2-4)

式(2-1)～(2-4)中，max{S_i,S_j}表示从S_i和S_j中取较大值，max{S_i,R-S_j×n}表示从S_i和R-S_j×n中取较大值，max{M-S_i×m,S_j}表示从M-S_i×m和S_j中取较大值，max{M-S_i×m,R-S_j×n}表示从M-S_i×m和R-S_j×n中取较大值，m表示矩阵A的行数M除以矩阵A被按行划分的行数S_i的结果的向下取整值，n表示矩阵B的列数R除以矩阵B被按列划分的列数S_j的结果的向下取整值，S_i表示矩阵A被按行划分的行数，S_j表示矩阵B被按列划分的列数，M表示矩阵A的行数，N表示矩阵A的列数，R表示矩阵B的列数。

式(2-2)～(2-4)中，变量k₁和k₂的值由式(1-3)和(1-4)决定。

式(1-3)和(1-4)中各个参数的定义可参见式(2-1)～(2-4)，在此不再赘述。

参见式(2-1)～(2-4)可知，评估函数(2)的形式十分复杂，本实施例中仅考虑S_i＝S_j的情况。因此，根据变量k₁和k₂的值可能有下述情形：

情形1(k₁＝0,k₂＝0)：

在这种情况下，m＝M/S_i,n＝R/S_j，式(2)的评估函数可简化为式(2-5)。

式(2-5)中，T_1,3(S_i,S_j)表示子矩阵①乘子矩阵③的总计算时钟节拍数，M表示矩阵A的行数，N表示矩阵A的列数，R表示矩阵B的列数，S_i表示矩阵A被按行划分的行数。

情形2(k₁＝0,k₂＝1)：

在这种情况下，m＝M/S_i,式(2)的评估函数可简化为式(2-6)。

式(2-6)中，T_1,3(S_i,S_j)表示子矩阵①乘子矩阵③的总计算时钟节拍数，T_1,4(S_i,S_j)表示子矩阵①乘子矩阵④的总计算时钟节拍数，M表示矩阵A的行数，N表示矩阵A的列数，R表示矩阵B的列数，S_j表示矩阵B被按列划分的列数。

情形3(k₁＝1,k₂＝0)：

在这种情况下，与情形2类似，n＝R/S_j，式(2)的评估函数可简化为式(2-7)。

式(2-7)中，T_1,3(S_i,S_j)表示子矩阵①乘子矩阵③的总计算时钟节拍数，T_2,3(S_i,S_j)表示子矩阵②乘子矩阵③的总计算时钟节拍数，M表示矩阵A的行数，N表示矩阵A的列数，R表示矩阵B的列数，S_i表示矩阵A被按行划分的行数。

情形4(k₁＝1,k₂＝1)：

令S'_i＝M-S_i×m,S'_j＝R-S_j×n，为了更好地简化讨论并且不失一般性，假设S'_i>S'_j，式(2)的评估函数可简化为式(2-8)。

式(2-8)中，T_1,3(S_i,S_j)表示子矩阵①乘子矩阵③的总计算时钟节拍数，T_1,4(S_i,S_j)表示子矩阵①乘子矩阵④的总计算时钟节拍数，T_2,3(S_i,S_j)表示子矩阵②乘子矩阵③的总计算时钟节拍数，T_2,4(S_i,S_j)表示子矩阵②乘子矩阵④的总计算时钟节拍数。

对于情形1至3，当S_i的值小于或等于处理单元链中处理单元的个数时，评估函数可取得最小值，因此在这三种情况下，最优的矩阵A被按行划分的行数S_i应当小于或等于处理单元链中处理单元的个数P。对于情形4，无法直接求得评估函数的最小值，因此必须借助数学工具对其进行分析，求得最优分块。对于具体的矩阵规模，M，N，R的值是已知的，也就意味着评估函数是一元函数。本实施例中，具体是使用MATLAB数学工具对评估函数生成满足式(1)约束的可视化图像，以矩阵A被按行划分的行数S_i和矩阵B被按列划分的列数S_j相等为前提，可以发现评估函数的自变量(分块大小)只能取整数值，并且在整数点函数导数为零，说明评估函数存在很多局部最值，本实施例通过作图分析的方式确定最优的矩阵A被按行划分的行数S_i，通过可视化图像确定评估函数定义域内的全局最小值，得到分块大小集合内最优的矩阵A被按行划分的行数S_i。

以规模为M＝128、N＝1728、R＝169的矩阵乘为例，此时评估函数的表达式具体如式(2-9)所示，并且自变量S_i满足约束条件如式(2-10)所示：

式(2-9)中，S_i表示矩阵A被按行划分的行数，即最优的分块大小；表示对128/S_i的结果向下取整；max为取较大值函数。本实施例中，具体是使用MATLAB数学工具对评估函数生成满足式(1)约束的可视化图像如图4所示。从图4可知，最优的矩阵A被按行划分的行数S_i在取值43时，评估函数达到最小值。矩阵A被按行划分的行数S_i的选取对矩阵加速器的性能发挥影响很大，只有唯一的一个最优的矩阵A被按行划分的行数S_i能使矩阵加速器发挥最高计算效率。

本实施例中，步骤3)的详细步骤包括：

3.1)将矩阵A以规模为S_i*N的子块为单位按行分块，得到包含整数个规模为S_i*N的子块的子矩阵①，如果仍有剩余的不规则的子块，则将剩余的不规则的子块作为子矩阵②；将矩阵B以规模为N*S_i的子块为单位按行分块，得到包含整数个规模为N*S_i的子块的子矩阵③，如果仍有剩余的不规则的子块，则将剩余的不规则的子块作为子矩阵④；将矩阵A和矩阵B的矩阵乘运算视为得到的所有子矩阵的两两相乘，根据得到的子矩阵数量建立使用指针相连接的DMA描述符链表，使得每一对相乘的子矩阵对应一个DMA描述符链表；

3.2)选择一个子块乘运算作为当前子块乘运算；

3.3)获取当前子块乘运算在矩阵A中对应子块数据的首地址、传输长度及传输步长，获取当前子块乘运算在矩阵B中对应子块数据的首地址、传输长度及传输步长，将当前子块乘运算在矩阵A中对应子块数据的首地址、传输长度及传输步长和当前子块乘运算在矩阵B中对应子块数据的首地址、传输长度及传输步长封装生成一个DMA描述符，根据当前子块乘运算所属的一对相乘的子矩阵确定对应的DMA描述符链表，并将该DMA描述符写入对应的DMA描述符链表中；

3.4)判断是否所有子块乘运算已经完成处理，如果尚未完成所有子块乘运算的处理，则选择下一个子块乘运算作为当前子块乘运算，跳转执行步骤3.3)；否则，如果已经完成所有子块乘运算的处理则跳转执行步骤4)。

假设矩阵A和矩阵B是两个规模为4*4的矩阵，当最优的子块大小S_i＝2时，将矩阵A分块为包含整数个规模为S_i*N的子块的子矩阵①、没有剩余的不规则的子矩阵②，将矩阵B分块为包含整数个规模为N*S_i的子块的子矩阵③、没有剩余的不规则的子矩阵④，因此只存在一对相乘的子矩阵子矩阵①乘子矩阵③，此时矩阵共被划分为四个子块，矩阵乘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分块为包含整数个规模为S_i*N的子块的子矩阵①和剩余的不规则的子矩阵②、将矩阵B分块为包含整数个规模为N*S_i的子块的子矩阵③和剩余的不规则的子矩阵④，因此共包含子矩阵①乘子矩阵③、子矩阵①乘子矩阵④、子矩阵②乘子矩阵③、子矩阵②乘子矩阵④共四对相乘的子矩阵，此时则需要生成四个DMA描述符链表，相邻DMA描述符链表之间使用指针相连接的DMA描述符链表。

本实施例中，步骤4)的详细步骤包括：

4.1)CPU将DMA描述符链表中第一个DMA描述符的首地址配置给矩阵乘加速器的DMA，矩阵乘加速器的DMA开始根据首地址读取第一个DMA描述符作为当前描述符；

4.2)矩阵乘加速器的DMA解析当前描述符，得到当前描述符对应子块乘运算所需的一对子块在外存中的地址及传输长度，并根据预设的传输步长和当前描述符中携带的一对子块的地址及传输长度，以分时的方式交换读取当前描述符对应子块乘运算所需数据，并将输入数据存入FIFO缓存；

4.3)基于FIFO缓存中当前描述符对应子块乘运算所需数据，通过矩阵乘加速器中至少一条处理单元链中的前S_i个处理单元进行当前描述符对应子块乘运算；

4.4)将当前描述符对应子块乘运算的结果通过矩阵乘加速器的DMA分别写回主存；

4.5)判断DMA描述符链表中的所有DMA描述符是否已经处理完毕，如果尚未处理完毕，则从DMA描述符链表中选择下一个DMA描述符作为当前描述符，跳转执行步骤4.2)；如果已经处理完毕，则结束并退出。

本实施例中，步骤4.3)的详细步骤包括：

4.3.1)在矩阵乘加速器中确定一条用于当前描述符对应子块乘运算的处理单元链，处理单元链至少包含S_i个处理单元；矩阵乘加速器中设有用于控制各个处理单元工作状态的状态机，状态机包括预取阶段、预取及计算阶段、提交阶段共三个状态；将变量k初始化为0，控制状态机进入预取状态；

4.3.2)在预取阶段，针对当前描述符对应子块乘运算，矩阵乘加速器的DMA将属于矩阵A的子块中的第k列数据附带上S_i个处理单元的编号信息后发出给处理单元链，处理单元链的前S_i个处理单元根据编号信息预取编号对应的数据，控制状态机进入预取及计算阶段；

4.3.3)在预取及计算阶段，矩阵乘加速器的DMA将属于矩阵B的子块中第k行数据以数据流的形式依次通过处理单元链的前S_i个处理单元，且前S_i个处理单元每一拍将接收到的数据和属于矩阵A的子块中的第k列数据部分进行乘法运算，并将乘法运算结果和片上存储中存储的上一次累加结果进行累加，将累加结果写回片上存储作为下一次累加的输入；同时，处理单元链的前S_i个处理单元分别根据编号信息预取属于矩阵A的子块中的第k+1列数据；

4.3.4)判断变量k的值是否等于N减1，如果变量k的值等于N减1，则控制状态机进入提交阶段，跳转执行步骤4.3.5)；否则，将变量k加1，控制状态机进入预取及计算阶段，跳转执行步骤4.3.3)；

4.3.5)在提交阶段下，处理单元链的前S_i个处理单元分别将最后一次累加的结果递给矩阵乘加速器的DMA，最终由矩阵乘加速器的DMA将结果写入主存中的指定区域。

如图6所示，本实施例在步骤3)中，CPU根据矩阵信息生成各个子块乘运算对应的DMA描述符(Buffer Descriptor,BD)，BD组织成四条DMA描述符链表存于主存中；DMA根据当前BD的信息就能够知道下一个BD在内存的位置，因此DMA只需知道第一个BD的首地址即可完成对所有BD的读取，每个BD包含了对应子块乘法的操作数的读取信息，DMA正是根据这些信息读取相应的数据。在整个读数过程中，CPU无需参与，只需等待最终结果写回即可，这极大减少了主机和DMA的通信开销。在进入步骤4)后，首先CPU将BD0(第一条DMA描述符链表的第一个DMA描述符)的首地址配置给DMA并启动DMA。DMA在工作状态下，通过数据流的形式，自动读取DMA描述符链表中的DMA描述符至FIFO中，当读到第四个DMA描述符链表BD4的最后一个DMA描述符时，DMA读取结束，等待矩阵加速器将结果写回主存。

如图7所示，本实施例中包含矩阵加速器的加速器系统基于xilinx公司的ZynqXC7Z045芯片上实现，主要包括处理机系统、片外存储和矩阵乘加速器，其中处理机系统包括CPU、AXI总线模块和存储控制器，CPU具体采用ARM双核A9处理器，CPU可以通过AXI总线模块经由存储控制器访问主存。处理机系统为矩阵加速器提供了AXI总线高速访存接口，矩阵加速器通过该接口也能访问主存，并与CPU共享主存空间。矩阵加速器由DMA和处理单元(Processing Element，PE)链(PE0～PEn)组成，每一个DMA和一个处理单元链相连，且由链首的PE0直接和DMA交互，DMA直接与高速访存接口相连接，可直接从主存读取加速器所需的数据。处理单元链是由若干个处理单元(PE)组成的链式结构，每个处理单元的结构完全一致，数据在相邻的处理单元间进行传递。初始化的时候处理数据存储在片外存储中，计算过程中的中间数据存在片上存储中，最终的结果又写回片外存储。具体过程是第一个处理单元PE0从DMA读取数据，依次向后一个PE传递数据，每个PE保存中间结果，等到全部计算完毕后，计算结果以相反方向向前一个PE传递，并最终由DMA将结果写回主存。矩阵加速器可扩展为若干条处理单元链，每条处理单元链都有其对应的DMA。由于高性能接口之间相对独立，都能够提供稳定的带宽，因此在DMA协助下，处理单元链之间可以高效并行地完成计算任务。

DMA主要用以管理片外存储于加速器间的数据传输。它的主要特点是计算大规模矩阵时往往只需要CPU配置一次，配置信息足够完整并存于片外存储中，DMA自动读取配置信息后就能够完成所有的数据传输工作。参见图6，DMA的工作步骤主要包括：步骤1，CPU根据矩阵分块后的矩阵的数据首地址、传输长度等信息生成DMA描述符，并将其组织成为链表结构，存入片外存储。一个矩阵子块乘法对应一个描述符，描述符之间用指针相连接。步骤2，CPU将描述符链表的在外存的首地址配置给DMA。DMA根据该地址自动读取DMA描述符，然后解析出描述符的配置信息，得到矩阵输入数据在外存的地址，并根据传输步长，传输长度等信息，以分时方式交换读取矩阵的输入数据。获得输入数据后，DMA将其存入相应的FIFO缓存，为加速器提供操作数据。步骤3，DMA读取完一个子块乘法所需的操作数之后，可以根据下个描述符的地址(当前描述符的信息)继续读取下一个子块乘法的操作数，在此过程中，DMA还负责将上一个子块乘法的运算结果写回外存。步骤4，重复步骤3直到DMA读取到最后一个DMA描述符。

如图8所示，本实施例中每个处理单元由计算模块、计算数据传送模块、数据存储模块和阶段同步控制模块组成。

计算单元采用Xilinx公司提供的浮点运算IP核进行实现。计算单元包括乘法器和加法器，其中加法器的流水线级别为11。

计算数据传送模块由FIFO_A、FIFO_B、FIFO_C共三个FIFO队列组成，FIFO_A和FIFO_B负责将输入数据传递到下一个PE。FIFO_C负责将结果传递到前一个PE中。

数据存储模块主要包括一个双端口BRAM(MEM_C)和地址生成器。当乘法器产生第一个结果时，地址生成器产生读信号和读地址，从MEM_C读取加法器的另一个操作数，也就是上一阶段的临时数据送往加法器；当加法器的第一个结果产生时，地址生成器产生写信号和写地址，结果被再次存入MEM_C。

阶段同步控制模块内部主要实现了两个计数器，当矩阵A子块的新的一列和矩阵B子块新的一行进入PE时，两个计数器分别开始计数。每个计数器都有相应的阈值，一旦其中一个计数器达到阈值，阶段同步控制模块就会使当前PE停止从前一个PE读取相对应的数据。这样一来达到阈值的计数器就会停止计数，一直到另一个计数器也达到其阈值，两个计数器同时被清零，重新开始为下一阶段的数据计数。通过阶段同步保证了矩阵A的某一列的第一个数据总是和矩阵B对应行的第一个数据同时进入每个PE。每个PE都有一个ID(PID)，矩阵A的每个数据在进入矩阵链之前被加入了编号信息用以表征这个数据归属于哪个处理单元。阶段同步控制模块，用以解决支持可变分块引起的数据不同步问题；此外，每个处理单元还包括用以控制与相邻处理单元间的数据交互的控制逻辑。

需要说明的是，矩阵加速器中处理单元的结构实现并不局限于如图8所示的特定结构，毫无疑问，本领域技术人员也可以根据需要采用其他结构的处理单元来实现矩阵加速器中的子块乘法运算，故在此不再赘述。

根据图2的原理可知，编号信息依次为0,1，···，S_i-1.当S_i<P时，那些满足PID>＝S_i的PE不允许启动工作。我们通过对数据编号和PID进行对比，使得当数据编号小于当前PID时，数据才允许被写入FIFO_A并传递到下一个PE。这样一来，数据就只传递到前S_i个PE中，我们的矩阵乘法器就能支持可变大小的分块。

本实施例中，每个处理单元内部实现了一个状态机用来控制数据流。处理单元工作时主要有三个阶段，我们以处理单元链中的第一个处理单元(PE0)为例：第一，预取阶段。PE0从DMA的FIFO中读取矩阵A的第一列数据，并通过FIFO_A传递给相邻的PE。数据从FIFO_A出来后在高位加入了编号信息，每个PE根据自己的ID号预取各自的数据。第二，计算-预取阶段。PE0开始同时读取矩阵A的第k+1列以及矩阵B的第k行数据并通过FIFO_A和FIFO_B传递给相邻的PE；矩阵B的数据与预取的A数据相乘，同时更新A的另一个预取缓冲(双缓冲，TA0，TA1)。经过一定延迟后对乘法结果进行加法运算，其中另一个操作数来自MEM_C，加法计算完毕后，又将结果存入MEM_C。第三阶段，提交阶段。每个PE将最后一次累加结果存入FIFO_C，同时MEM_C内容清零。存入完毕后，每个PE从后一个PE读取结果存入FIFO_C，由此结果数据依次前移至DMA的FIFO_C中。

本实施例以典型卷积神经网络(CNN)中如表1所示的五个非均匀矩阵规模Conv1～Conv5为例，应用本实施例以及经典分块算法分别对表1中的五个规模的矩阵乘进行加速对比，表2是本实施例中每个卷积层的矩阵规模通过优化分块策略算出的最优的矩阵A被按行划分的行数S_i。

表1：五个非均匀矩阵规模实例表。

卷积层	M	N	R
				Conv1	96	363	3025
Conv2	128	1200	729
				Conv3	384	2304	169
Conv4	192	1728	169
				Conv5	128	1728	169

表2：五个非均匀矩阵规模实例的最优的矩阵A被按行划分的行数S_i数据表。

最终，得到的性能对比图如图9所示。参见图9可知，当矩阵规模不均匀程度比较高时(矩阵的行列大小相差较大)，本实施例的矩阵乘甲酸方法能达到很高的性能提升，例如对于Conv5,本实施例的矩阵乘甲酸方法相对于经典分块算法有大约12％的性能提升。

以上所述仅是本发明的优选实施方式，本发明的保护范围并不仅局限于上述实施例，凡属于本发明思路下的技术方案均属于本发明的保护范围。应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理前提下的若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。

Claims (4)

1.一种支持可变分块的矩阵乘加速方法，其特征在于步骤包括：

1)输入矩阵乘运算所需的M*N的矩阵A和N*R的矩阵B；

2)根据矩阵A和矩阵B的规模确定矩阵A被按行划分的行数S_i，将矩阵A以规模为S_i*N的子块为单位进行按行分块，将矩阵B以规模为N*S_i的子块为单位进行按列分块，使得矩阵乘运算等同为多个子块乘运算；

3)为每一个子块乘运算所需数据生成一个DMA描述符，将所有子块乘运算的DMA描述符构建DMA描述符链表并存入主存；

4)针对每一个子块乘运算，通过矩阵乘加速器的DMA从主存的DMA描述符链表读取子块乘运算所需数据，然后通过矩阵乘加速器中至少一条处理单元链中的前S_i个处理单元进行子块乘运算，并将各个子块乘运算的结果通过矩阵乘加速器的DMA分别写回主存；

所述步骤2)中根据矩阵A和矩阵B的规模确定矩阵A被按行划分的行数S_i的详细步骤包括：

2.1)输入矩阵A和矩阵B的矩阵规模，所述矩阵规模包括矩阵A的行数M、矩阵A的列数N、矩阵B的列数R三者的值；

2.2)根据矩阵规模获取满足式(1)所示约束的分块大小集合；

<mrow> <mfenced open = "{" close = ""> <mtable> <mtr> <mtd> <mrow> <mo>(</mo> <msub> <mi>S</mi> <mi>i</mi> </msub> <mo>+</mo> <msub> <mi>S</mi> <mi>j</mi> </msub> <mo>)</mo> <mo>&amp;times;</mo> <mi>N</mi> <mo>/</mo> <mi>B</mi> <mi>W</mi> <mo>&amp;le;</mo> <mo>(</mo> <mn>2</mn> <mo>&amp;times;</mo> <msub> <mi>S</mi> <mi>i</mi> </msub> <mo>+</mo> <mi>m</mi> <mi>a</mi> <mi>x</mi> <mo>{</mo> <msub> <mi>S</mi> <mi>i</mi> </msub> <mo>,</mo> <msub> <mi>S</mi> <mi>j</mi> </msub> <mo>}</mo> <mo>&amp;times;</mo> <mi>N</mi> <mo>)</mo> <mo>/</mo> <mi>F</mi> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <msub> <mi>S</mi> <mi>i</mi> </msub> <mo>&lt;</mo> <mo>=</mo> <mi>P</mi> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <mi>m</mi> <mi>a</mi> <mi>x</mi> <mo>{</mo> <msub> <mi>S</mi> <mi>i</mi> </msub> <mo>,</mo> <msub> <mi>S</mi> <mi>j</mi> </msub> <mo>}</mo> <mo>&gt;</mo> <msub> <mi>Stage</mi> <mrow> <mi>a</mi> <mi>d</mi> <mi>d</mi> </mrow> </msub> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <mi>m</mi> <mi>a</mi> <mi>x</mi> <mo>{</mo> <mi>M</mi> <mo>-</mo> <msub> <mi>S</mi> <mi>i</mi> </msub> <mo>&amp;times;</mo> <mi>m</mi> <mo>,</mo> <mi>R</mi> <mo>-</mo> <msub> <mi>S</mi> <mi>j</mi> </msub> <mo>&amp;times;</mo> <mi>n</mi> <mo>}</mo> <mo>&gt;</mo> <msub> <mi>Stage</mi> <mrow> <mi>a</mi> <mi>d</mi> <mi>d</mi> </mrow> </msub> </mrow> </mtd> </mtr> </mtable> </mfenced> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>1</mn> <mo>)</mo> </mrow> </mrow>

式(1)中，S_i表示矩阵A被按行划分的行数，S_j表示矩阵B被按列划分的列数，M表示矩阵A的行数，N表示矩阵A的列数，BW表示矩阵乘加速器中单条处理单元链的访存带宽，max{S_i,S_j}表示从S_i和S_j中取较大值，F表示访存频率；P表示矩阵乘加速器中单条处理单元链的处理单元数量，Stage_add表示矩阵乘加速器中的加法器流水级数，max{M-S_i×m,R-S_j×n}表示从M-S_i×m和R-S_j×n中取较大值，R表示矩阵B的列数，m表示矩阵A的行数M除以矩阵A被按行划分的行数S_i的结果的向下取整值，n表示矩阵B的列数R除以矩阵B被按列划分的列数S_j的结果的向下取整值；

2.3)将矩阵A视为包含整数个规模为S_i*N的子块的子矩阵①和剩余的不规则的子矩阵②、将矩阵B视为包含整数个规模为N*S_i的子块的子矩阵③和剩余的不规则的子矩阵④，将矩阵A和矩阵B的矩阵乘运算视为子矩阵①～子矩阵④四者中的两两相乘，建立式(2)所示的评估函数；

f(S_i,S_j)＝T_1,3(S_i,S_j)+T_1,4(S_i,S_j)+T_2,3(S_i,S_j)+T_2,4(S_i,S_j) (2)

式(2)中，f(S_i,S_j)表示评估函数的值，S_i表示矩阵A被按行划分的行数，S_j表示矩阵B被按列划分的列数，T_1,3(S_i,S_j)表示子矩阵①乘子矩阵③的计算时钟节拍数，T_1,4(S_i,S_j)表示子矩阵①乘子矩阵④的计算时钟节拍数，T_2,3(S_i,S_j)表示子矩阵②乘子矩阵③的计算时钟节拍数，T_2,4(S_i,S_j)表示子矩阵②乘子矩阵④的计算时钟节拍数；

2.4)使用MATLAB数学工具对式(2)所示的评估函数生成满足式(1)约束的可视化图像，以矩阵A被按行划分的行数S_i和矩阵B被按列划分的列数S_j相等为前提，通过所述可视化图像确定评估函数定义域内的全局最小值，得到分块大小集合内最优的矩阵A被按行划分的行数S_i。

2.根据权利要求1所述的支持可变分块的矩阵乘加速方法，其特征在于，所述步骤3)的详细步骤包括：

3.1)将矩阵A以规模为S_i*N的子块为单位按行分块，得到包含整数个规模为S_i*N的子块的子矩阵①，如果仍有剩余的不规则的子块，则将剩余的不规则的子块作为子矩阵②；将矩阵B以规模为N*S_i的子块为单位按行分块，得到包含整数个规模为N*S_i的子块的子矩阵③，如果仍有剩余的不规则的子块，则将剩余的不规则的子块作为子矩阵④；将矩阵A和矩阵B的矩阵乘运算视为得到的所有子矩阵的两两相乘，根据得到的子矩阵数量建立使用指针相连接的DMA描述符链表，使得每一对相乘的子矩阵对应一个DMA描述符链表；

3.2)选择一个子块乘运算作为当前子块乘运算；

3.3)获取当前子块乘运算在矩阵A中对应子块数据的首地址、传输长度及传输步长，获取当前子块乘运算在矩阵B中对应子块数据的首地址、传输长度及传输步长，将当前子块乘运算在矩阵A中对应子块数据的首地址、传输长度及传输步长和当前子块乘运算在矩阵B中对应子块数据的首地址、传输长度及传输步长封装生成一个DMA描述符，根据当前子块乘运算所属的一对相乘的子矩阵确定对应的DMA描述符链表，并将该DMA描述符写入对应的DMA描述符链表中；

3.4)判断是否所有子块乘运算已经完成处理，如果尚未完成所有子块乘运算的处理，则选择下一个子块乘运算作为当前子块乘运算，跳转执行步骤3.3)；否则，如果已经完成所有子块乘运算的处理则跳转执行步骤4)。

3.根据权利要求2所述的支持可变分块的矩阵乘加速方法，其特征在于，所述步骤4)的详细步骤包括：

4.1)CPU将DMA描述符链表中第一个DMA描述符的首地址配置给矩阵乘加速器的DMA，矩阵乘加速器的DMA开始根据首地址读取第一个DMA描述符作为当前描述符；

4.2)矩阵乘加速器的DMA解析当前描述符，得到当前描述符对应子块乘运算所需的一对子块在外存中的地址及传输长度，并根据预设的传输步长和当前描述符中携带的一对子块的地址及传输长度，以分时的方式交换读取当前描述符对应子块乘运算所需数据，并将输入数据存入FIFO缓存；

4.3)基于FIFO缓存中当前描述符对应子块乘运算所需数据，通过矩阵乘加速器中至少一条处理单元链中的前S_i个处理单元进行当前描述符对应子块乘运算；

4.4)将当前描述符对应子块乘运算的结果通过矩阵乘加速器的DMA分别写回主存；

4.5)判断DMA描述符链表中的所有DMA描述符是否已经处理完毕，如果尚未处理完毕，则从DMA描述符链表中选择下一个DMA描述符作为当前描述符，跳转执行步骤4.2)；如果已经处理完毕，则结束并退出。

4.根据权利要求3所述的支持可变分块的矩阵乘加速方法，其特征在于，所述步骤4.3)的详细步骤包括：

4.3.1)在矩阵乘加速器中确定一条用于当前描述符对应子块乘运算的处理单元链，所述处理单元链至少包含S_i个处理单元；所述矩阵乘加速器中设有用于控制各个处理单元工作状态的状态机，所述状态机包括预取阶段、预取及计算阶段、提交阶段共三个状态；将变量k初始化为0，控制状态机进入预取状态；

4.3.2)在预取阶段，针对当前描述符对应子块乘运算，矩阵乘加速器的DMA将属于矩阵A的子块中的第k列数据附带上S_i个处理单元的编号信息后发出给处理单元链，处理单元链的前S_i个处理单元根据编号信息预取编号对应的数据，控制状态机进入预取及计算阶段；

4.3.3)在预取及计算阶段，矩阵乘加速器的DMA将属于矩阵B的子块中第k行数据以数据流的形式依次通过处理单元链的前S_i个处理单元，且前S_i个处理单元每一拍将接收到的数据和属于矩阵A的子块中的第k列数据部分进行乘法运算，并将乘法运算结果和片上存储中存储的上一次累加结果进行累加，将累加结果写回片上存储作为下一次累加的输入；同时，处理单元链的前S_i个处理单元分别根据编号信息预取属于矩阵A的子块中的第k+1列数据；

4.3.4)判断变量k的值是否等于N减1，如果变量k的值等于N减1，则控制状态机进入提交阶段，跳转执行步骤4.3.5)；否则，将变量k加1，控制状态机进入预取及计算阶段，跳转执行步骤4.3.3)；

4.3.5)在提交阶段下，处理单元链的前S_i个处理单元分别将最后一次累加的结果递给矩阵乘加速器的DMA，最终由矩阵乘加速器的DMA将结果写入主存中的指定区域。