CN104346318A

CN104346318A - 面向通用多核dsp的矩阵乘加速方法

Info

Publication number: CN104346318A
Application number: CN201410541958.1A
Authority: CN
Inventors: 迟利华; 刘杰; 甘新标; 晏益慧; 徐涵; 胡庆丰; 蒋杰; 李胜国; 王庆林; 皇甫永硕; 崔显涛; 周陈
Original assignee: National University of Defense Technology
Current assignee: National University of Defense Technology
Priority date: 2014-10-15
Filing date: 2014-10-15
Publication date: 2015-02-11
Anticipated expiration: 2034-10-15
Also published as: CN104346318B

Abstract

本发明公开了一种面向通用多核DSP的矩阵乘加速方法，目的是提高矩阵乘的运算速度和最大化通用多核DSP的计算效率。技术方案是先对进行DSP配置并初始化：然后对矩阵A和矩阵B进行划分，依据VPU的拓扑结构mg×ng将原矩阵乘运算转换为分块矩阵乘运算；接着每个VPU同时并行执行数据迁移操作；最后依据数据分布原理，将各VPU中AM的结果归并组成结果矩阵C＝A×B的计算结果。采用本发明可以提高通用多核DSP结构的矩阵乘运算速度和通用多核DSP系统的计算资源利用率。

Description

面向通用多核DSP的矩阵乘加速方法

技术领域

本发明涉及矩阵乘加速方法，尤指面向通用多核DSP结构的矩阵乘加速方法。

背景技术

随着通用DSP计算性能的不断攀升和通用DSP的广泛应用，通用多核DSP必将成为高性能计算重要的发展方向。矩阵乘是数值计算中最常用的一类操作，很多应用中都包含矩阵乘的计算过程，面向通用多核DSP设计高效的矩阵乘方法可以有效提高应用的计算速度和提升通用多核DSP的计算效率，以达到通用多核DSP的设计目标。

矩阵乘是将被乘矩阵A的一行和乘数矩阵B的一列相乘得到结果矩阵C中的一个元素。面向多核系统的矩阵乘法通常需要将矩阵乘计算过程合理分布在各个处理器计算核之间并行完成计算过程，以提高矩阵乘的运算速度和最大化异构系统的计算效率。

由于处理器设计目标和指令集结构的不同，传统的面向通用处理器的矩阵乘实现技术很难满足面向特定应用而设计的处理器DSP(Digital Signal Processing,数字信号处理)的性能要求，因此，必须面向DSP体系结构定制高效的矩阵乘，以提高矩阵乘的运算速度，最大限度满足处理器设计目标。随着通用多核DSP计算性能的不断攀升和广泛应用，面向通用多核DSP定制高效矩阵乘是DSP在高性能计算领域发展的关键技术。

高效的矩阵乘与处理器体系结构紧密相关。因此，面向通用多核DSP体系结构，如图1所示，如何面向多核DSP结构对矩阵乘进行加速是本领域技术人员亟需解决的技术难题。

发明内容

本发明要解决的技术问题在于：面向通用多核DSP体系结构提出一种面向通用多核DSP的矩阵乘加速方法，以提高矩阵乘的运算速度和最大化通用多核DSP的计算效率。

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

第一步、DSP配置与初始化。

1.1定义矩阵A的维度为M×K，矩阵B的维度为K×N，则A和B相乘的结果矩阵C的维度为M×N，M,K,N均为正整数；A的元素为a_pq(0≤p≤M-1,0≤q≤K-1)，B的元素为b_st(0≤s≤K-1,0≤t≤N-1)；

1.2查询体系结构手册获取DSP中DSP计算核，即虚拟处理单元VPU(VirtualProcessing Unit)的拓扑结构mg×ng，即DSP中有mg×ng个VPU，物理分布为mg行×ng列，依次编号为(0,0),(0,1),…(0,ng-1)，(1,0),(1,1),…(1,ng-1)，……(mg-1,0),(mg-1,1),…(mg-1,ng-1)；

1.3查询体系结构手册获取每个VPU拥有的浮点乘累加功能部件VPE(VirtualProcessing Element)的数目m_e；

1.4利用系统软件提供的初始化函数(如init)完成DSP初始化；

第二步、对矩阵A和矩阵B进行划分，依据VPU的拓扑结构mg×ng将原矩阵乘运算转换为分块矩阵乘运算，具体方法如下：

2.1将M*K的矩阵A采用矩阵划分方法划分为mg×ng个A块矩阵，A块矩阵表示为A_ij(0≤i≤mg-1,0≤j≤ng-1)，每个块矩阵的维度为m×k，其中1≤m≤M，1≤k≤K，mg*m+mod(M,mg)＝M，其中，mod(M,mg)表示对M用mg求余，ng*k+mod(K,ng)＝K，其中，mod(K,ng)表示对K用ng求余；边缘A块矩阵的维度需要特殊处理。具体的矩阵划分方法如下：

2.1.1令A块矩阵行变量i＝0；

2.1.2令A块矩阵列变量j＝0；

2.1.3令A行坐标变量s＝0；

2.1.4令A列坐标变量e＝0；

2.1.5令A块矩阵元素行变量m'＝0；

2.1.6令A块矩阵元素列变量n'＝0；

2.1.7s＝m*i+m'；

2.1.8e＝n*j+n'；

2.1.9如果s≤M-1，转2.1.10，否则，矩阵A划分完毕，转2.2；

2.1.10如果e≤K-1，转2.1.11，否则，矩阵A划分完毕，转2.2；

2.1.11选取矩阵A的元素a_se组成A块矩阵Ai_j的元素a'_m'n'；

2.1.12n'＝n'+1；

2.1.13如果n'≤k-1，转2.1.7，否则，转2.1.14；

2.1.14m'＝m'+1；

2.1.15如果m'≤m-1，转2.1.6，否则，转2.1.16；

2.1.16j＝j+1；

2.1.17如果j≤ng-1，转2.1.3，否则，转2.1.18；

2.1.18i＝i+1；

2.1.19如果i≤mg-1，转2.1.2，否则，矩阵A划分完毕，转2.2；

2.2将K*N的矩阵B采用矩阵划分方法划分为mg×ng个B块矩阵，B块矩阵表示为B_pq，0≤p≤mg-1,0≤q≤ng-1，每个B块矩阵的维度为k×n，其中1≤k≤K,1≤n≤N，mg*k+mod(K,mg)＝K，其中，mod(K,mg)表示求对K用mg求余，ng*n+mod(K,ng)＝N，其中，mod(K,ng)表示对K用ng求余；边缘B块矩阵的维度需要特殊处理，具体的矩阵划分方法如下：

2.2.1令B块矩阵行变量p＝0；

2.2.2令B块矩阵列变量q＝0；

2.2.3令B行坐标变量s'＝0；

2.2.4令B列坐标变量e'＝0；

2.2.5令B块矩阵元素行变量k'＝0；

2.2.6令B块矩阵元素列变量n″＝0；

2.2.7s'＝k*p+k'；

2.2.8e'＝n*q+n″；

2.2.9如果s'≤K-1，转2.2.10，否则，矩阵B划分完毕，转2.3；

2.2.10如果e'≤N-1，转2.2.11，否则，矩阵B划分完毕，转2.3；

2.2.11选取矩阵B的元素b_s'e'组成块矩阵B_pq的元素b'_k'n″；

2.2.12n″＝n″+1；

2.2.13如果n″≤n-1，转2.2.7，否则，转2.2.14；

2.2.14k'＝k'+1；

2.2.15如果k'≤k-1，转2.2.6，否则，转2.2.16；

2.2.16q＝q+1；

2.2.17如果j2≤ng-1，转2.2.3，否则，转2.2.18；

2.2.18p＝p+1；

2.2.19如果i2≤mg-1，转2.2.2，否则，矩阵B划分完毕，转2.3；

2.3初始化结果矩阵C＝0，即令C中所有元素都为0，即元素C(m₀,n₀)＝0，1≤m0≤M,1≤n0≤N。

第三步、由于每个VPU是相互独立的处理单元，每个VPU同时并行执行如下数据迁移操作：

3.1令A块矩阵行变量i'＝0；

3.2令A块矩阵列变量j'＝0；

3.3核(i',j')从DSP存储空间读A块矩阵A_i'j'到阵列存储器(Array Memory)，A_i'j'占用的存储空间记为

3.4令B块矩阵列变量k″＝0；

3.5核(j',k″)从DSP存储空间读B块矩阵B_j'k″到标量数据缓冲区L1D；

3.6由于每个VPU是相互独立的计算单元，每个VPU并行执行如下向量乘法操作：

3.7获取A块矩阵A_i'j'的维度m'0×k0，1≤m'0≤m,1≤k0≤k；

3.8令列长度变量j″＝0；

3.9初始化向量循环变量v＝0；

3.10定义向量循环次数表示下取整；

3.11初始化广播变量r＝0；

3.12从A块矩阵A_i'j'的第j″列的第v*m_e个元素开始连续读取m_e个元素，组成m_e维的向量

V_{a} = (a_{(v m_{e} + 0) j^{''}}, a_{({vm}_{e} + 1) j^{''}}, \cdot \cdot \cdot, a_{({vm}_{e} + (m_{e} - 1)) j^{''}});

3.13取矩阵B_j'k″第j″行中的第r个元素广播成一个m_e维的向量

3.14向量V_a与V_b对应元素相乘，得到矩阵C第j″列连续m_e个元素的部分乘积，具体相乘规则如下：

3.15将3.14得到的矩阵C第j″列连续m_e个元素的部分乘积传回至阵列存储器，这m_e个元素占据的存储空间记为

3.16如果j″＞0，转3.17，否则，转3.18；

3.17将3.15传回的第j″列连续m_e个元素的部分乘积与AM中存储的j″-1列连续m_e个元素的部分乘积对应相加，具体相加方法如下：

3.17.1令i3＝vm_e+0；

3.17.2C_i'k″(i3,j″)＝C_i'k″(i3,j″)+C_i'k″(i3,j″-1)，C_i'k″(i3,j″)表示结果矩阵的块矩阵C_i'k″的元素；

3.17.3i3＝i3+1；

3.17.4如果i3≤vm_e+(m_e-1)，转3.17.2，否则，转3.17.5；

3.17.5将3.17.2累加结果传回至阵列存储器，占用的存储空间记为

3.18广播变量r＝r+1；

3.19如果r＜k0，转3.12，否则，转3.21；

3.20向量循环变量v＝v+1；

3.21如果v＜Vz，转3.11，否则，转3.23；

3.22定义向量循环余数vr＝m0-v*m_e；

3.23从A块矩阵A_i'j'的第j″列的第v*m_e个元素开始连续读取vr个元素构成一个m_e维的向量的前vr个分量，后面m_e-vr个分量以0补充，组成向量

3.24如果v＝Vz，转3.13，否则，转3.26；

3.25j″＝j″+1；

3.26如果j″＜n0；转3.9，否则，转3.28；

3.27k″＝k″+1；

3.28如果k″≤ng-1，转3.5，否则，转3.30；

3.29j'＝j'+1；

3.30如果j'≤ng-1，转3.3，否则，转3.32；

3.31i'＝i'+1；

3.32如果i≤mg-1，转3.2，否则，转第四步；

第四步、依据数据分布原理，将各VPU中阵列存储器的结果归并组成结果矩阵C＝A×B的计算结果。具体方法如下：

4.1令矩阵C行变量u＝0；

4.2令矩阵C列变量v＝0；

4.3将块矩阵C_uv的计算结果传回至DSP全局存储器Global Memory；

4.4获取块矩阵C_uv的维度mc×nc，1≤mc≤m,1≤nc≤n；

4.5令行坐标变量i_c＝0；

4.6令列坐标变量j_c＝0；

4.7C(i_c,j_c)＝C(i_c,j_c)+C_uv(i_c,j_c)，C(i_c,j_c)表示矩阵C的元素，C_uv(i_c,j_c)表示块矩阵C_uv的元素；

4.8j_c＝j_c+1；

4.9如果j_c＜nc；转4.7，否则，转4.10；

4.10i_c＝i_c+1；

4.11如果i_c＜mc，转4.6，否则，转4.12；

4.12v＝v+1；

4.13如果v≤ng-1；转4.3，否则，转4.14；

4.14u＝u+1；

4.15如果u≤mg-1，转4.2，否则，归并完毕，转第五步；

第五步、结束。

采用本发明可以达到以下技术效果：

1.本发明第三步由每个VPU同时并行执行数据迁移和数据计算操作，提高了通用多核DSP结构的矩阵乘运算速度；

2.本发明第二步依据VPU的拓扑结构划分矩阵数据，提升了通用多核DSP系统的计算资源利用率。

附图说明

图1为通用多核DSP体系结构；

图2为本发明面向通用多核DSP的矩阵乘加速方法总体流程图。

具体实施方式

图1为通用多核DSP体系结构；

在图1中，每个单核DSP由SPU和VPU构成，SPU由L1I(Level 1 Instruction)Cache、L1D(Level 1 Data)Cache、SPE(Special Processing Unit)、标量寄存器和流控制器组成，L1I用于指令缓存；L1D用于数据缓存；SPE用于一些指令流的控制、对向量单元的配置以及主要的通信任务；VPU包括AM(Array Memory)、向量寄存器和多个可并发执行的VPE(Virtual Processing Unit)，AM主要用于数组缓存，多个VPE组成的向量SIMD单元主要用于数值运算加速。

图2为本发明面向通用多核DSP的矩阵乘加速方法总体流程图

本发明的步骤如下：

第一步、DSP配置与初始化。

第二步、对矩阵A和矩阵B进行划分，依据VPU的拓扑结构mg×ng将原矩阵乘运算转换为分块矩阵乘运算。

第三步、每个VPU同时并行执行数据迁移操作。

第四步、依据数据分布原理，将各VPU中AM的结果归并组成结果矩阵。

第五步、结束。

Claims

1.一种面向通用多核DSP的矩阵乘加速方法，其特征在于包括以下步骤：

第一步、DSP配置与初始化：

1.1 定义矩阵A的维度为M×K，矩阵B的维度为K×N，A和B相乘的结果矩阵C的维度为M×N，M,K,N均为正整数；A的元素为a_pq(0≤p≤M-1,0≤q≤K-1)，B的元素为b_st(0≤s≤K-1,0≤t≤N-1)；

1.2 查询体系结构手册获取DSP中DSP计算核，即虚拟处理单元VPU的拓扑结构mg×ng，即DSP中有mg×ng个VPU，物理分布为mg行×ng列，依次编号为(0，0)，(0，1)，…(0，ng-1)，(1，0)，(1，1)，…(1，ng-1)，……(mg-1，0)，(mg-1，1)，…(mg-1，ng-1)；

1.3 查询体系结构手册获取每个VPU拥有的浮点乘累加功能部件VPE的数目m_e；

1.4 利用系统软件提供的初始化函数完成DSP初始化；

第二步、对矩阵A和矩阵B进行划分，依据VPU的拓扑结构mg×ng将原矩阵乘运算转换为分块矩阵乘运算；

第三步、每个VPU同时并行执行如下数据迁移操作：

3.1 令A块矩阵行变量i'＝0；

3.2 令A块矩阵列变量j'＝0；

3.3 核(i',j')从DSP存储空间读A块矩阵A_i'j'到存储器AM，A_i'j'占用的存储空间记为

3.4 令B块矩阵列变量k″＝0；

3.5 核(j',k″)从DSP存储空间读B块矩阵B_j'k″到标量数据缓冲区L1D；

3.6 每个VPU并行执行如下向量乘法操作：

3.7 获取A块矩阵A_i'j'的维度m'0×k0，1≤m'0≤m,1≤k0≤k；

3.8 令列长度变量j″＝0；

3.9 初始化向量循环变量v＝0；

3.10 定义向量循环次数表示下取整；

3.11 初始化广播变量r＝0；

3.12 从A块矩阵A_i'j'的第j″列的第v*m_e个元素开始连续读取m_e个元素，组成m_e维的向量

V_{a} = (a_{({vm}_{e} + 0) j^{''}}, a_{({vm}_{e} + 1) j^{''}}, . . ., a_{({vm}_{e} + (m_{e} - 1)) j^{''}});

3.13 取矩阵B_j′k′第j″行中的第r个元素广播成一个m_e维的向量

3.14 向量V_a与V_b对应元素相乘，得到矩阵C第j″列连续m_e个元素的部分乘积，具体相乘规则如下：

3.15 将3.14得到的矩阵C第j″列连续m_e个元素的部分乘积传回至AM，这m_e个元素占据的存储空间记为

3.16 如果j″＞0，转3.17，否则，转3.18；

3.17 将3.15传回的第j″列连续m_e个元素的部分乘积与AM中存储的j″-1列连续m_e个元素的部分乘积对应相加，具体相加方法如下：

3.17.1 令i3＝vm_e+0；

3.17.2 C_i'k″(i3,j″)＝C_i'k″(i3,j″)+C_i'k″(i3,j″-1)，C_i'k″(i3,j″)表示结果矩阵的块矩阵C_i'k″的元素；

3.17.3 i3＝i3+1；

3.17.4 如果i3≤vm_e+(m_e-1)，转3.17.2，否则，转3.17.5；

3.17.5 将3.17.2累加结果传回至阵列存储器，占用的存储空间记为

3.18 广播变量r＝r+1；

3.19 如果r＜k0，转3.12，否则，转3.21；

3.20 向量循环变量v＝v+1；

3.21 如果v＜Vz，转3.11，否则，转3.23；

3.22 定义向量循环余数vr＝m0-v*m_e；

3.23 从A块矩阵A_i'j'的第j″列的第v*m_e个元素开始连续读取vr个元素构成一个m_e维的向量的前vr个分量，后面m_e-vr个分量以0补充，组成向量

3.24 如果v＝Vz，转3.13，否则，转3.26；

3.25 j″＝j″+1；

3.26 如果j″＜n0；转3.9，否则，转3.28；

3.27 k″＝k″+1；

3.28 如果k″≤ng-1，转3.5，否则，转3.30；

3.29 j'＝j'+1；

3.30 如果j'≤ng-1，转3.3，否则，转3.32；

3.31 i'＝i'+1；

3.32 如果i≤mg-1，转3.2，否则，转第四步；

第四步、依据数据分布原理，将各VPU中阵列存储器的结果归并组成结果矩阵C＝A×B的计算结果；

第五步、结束。

2.如权利要求1所述的面向通用多核DSP的矩阵乘加速方法，其特征在于对矩阵A和矩阵B进行划分的方法是：

2.1 将M*K的矩阵A采用矩阵划分方法划分为mg×ng个A块矩阵，A块矩阵表示为A_ij，0≤i≤mg-1,0≤j≤ng-1，每个块矩阵的维度为m×k，其中1≤m≤M,1k≤，≤mg*m+mod(M,mg)＝M，其中，mod(M,mg)表示对M用mg求余，ng*k+mod(K,ng)＝K，其中，mod(K,ng)表示对K用ng求余；边缘A块矩阵的维度需要特殊处理；具体的矩阵划分方法如下：

2.1.1 令A块矩阵行变量i＝0；

2.1.2 令A块矩阵列变量j＝0；

2.1.3 令A行坐标变量s＝0；

2.1.4 令A列坐标变量e＝0；

2.1.5 令A块矩阵元素行变量m'＝0；

2.1.6 令A块矩阵元素列变量n'＝0；

2.1.7 s＝m*i+m'；

2.1.8 e＝n*j+n'；

2.1.9 如果s≤M-1，转2.1.10，否则，矩阵A划分完毕，转2.2；

2.1.10 如果e≤K-1，转2.1.11，否则，矩阵A划分完毕，转2.2；

2.1.11 选取矩阵A的元素a_se组成A块矩阵A_ij的元素a'_m'n'；

2.1.12 n'＝n'+1；

2.1.13 如果n'≤k-1，转2.1.7，否则，转2.1.14；

2.1.14 m'＝m'+1；

2.1.15 如果m'≤m-1，转2.1.6，否则，转2.1.16；

2.1.16 j＝j+1；

2.1.17 如果j≤ng-1，转2.1.3，否则，转2.1.18；

2.1.18 i＝i+1；

2.1.19 如果i≤mg-1，转2.1.2，否则，矩阵A划分完毕，转2.2；

2.2 将K*N的矩阵B采用矩阵划分方法划分为mg×ng个B块矩阵，B块矩阵表示为B_pq，0≤p≤mg-1,0≤q≤ng-1，每个B块矩阵的维度为k×n，其中1≤k≤K,1≤n≤N，mg*k+mod(K,mg)＝K，其中，mod(K,mg)表示求对K用mg求余，ng*n+mod(K,ng)＝N，其中，mod(K,ng)表示对K用ng求余；边缘B块矩阵的维度需要特殊处理，具体的矩阵划分方法如下：

2.2.1 令B块矩阵行变量p＝0；

2.2.2 令B块矩阵列变量q＝0；

2.2.3 令B行坐标变量s'＝0；

2.2.4 令B列坐标变量e'＝0；

2.2.5 令B块矩阵元素行变量k'＝0；

2.2.6 令B块矩阵元素列变量n″＝0；

2.2.7 s'＝k*p+k'；

2.2.8 e'＝n*q+n″；

2.2.9 如果s'≤K-1，转2.2.10，否则，矩阵B划分完毕，转2.3；

2.2.10 如果e'≤N-1，转2.2.11，否则，矩阵B划分完毕，转2.3；

2.2.11 选取矩阵B的元素b_s'e'组成块矩阵B_pq的元素b'_k'n″；

2.2.12 n″＝n″+1；

2.2.13 如果n″≤n-1，转2.2.7，否则，转2.2.14；

2.2.14 k'＝k'+1；

2.2.15 如果k'≤k-1，转2.2.6，否则，转2.2.16；

2.2.16 q＝q+1；

2.2.17 如果j2≤ng-1，转2.2.3，否则，转2.2.18；

2.2.18 p＝p+1；

2.2.19 如果i2≤mg-1，转2.2.2，否则，矩阵B划分完毕，转2.3；

2.3 初始化结果矩阵C＝0，即令C中所有元素都为0，即元素C(m₀,n₀)＝0，1≤m0≤M,1≤n0≤N。

3.如权利要求1所述的面向通用多核DSP的矩阵乘加速方法，其特征在于第四步将各VPU中阵列存储器的结果归并组成结果矩阵C＝A×B的计算结果的具体方法如下：

3.1 令矩阵C行变量u＝0；

3.2 令矩阵C列变量v＝0；

3.3 将块矩阵C_uv的计算结果传回至DSP全局存储器Global Memory；

3.4 获取块矩阵C_uv的维度mc×nc，1≤mc≤m,1≤nc≤n；

3.5 令行坐标变量i_c＝0；

3.6 令列坐标变量j_c＝0；

3.7 C(i_c,j_c)＝C(i_c,j_c)+C_uv(i_c,j_c)，C(i_c,j_c)表示矩阵C的元素，C_uv(i_c,j_c)表示块矩阵C_uv的元素；

3.8 j_c＝j_c+1；

3.9 如果j_c＜nc；转4.7，否则，转4.10；

3.10 i_c＝i_c+1；

3.11 如果i_c＜mc，转4.6，否则，转4.12；

3.12 v＝v+1；

3.13 如果v≤ng-1；转4.3，否则，转4.14；

3.14 u＝u+1；

3.15 如果u≤mg-1，转4.2，否则，归并完毕。