面向cpu+dsp异构系统的矩阵乘加速方法
【专利摘要】本发明公开了一种面向CPU+DSP异构系统的矩阵乘加速方法，目的是面向CPU+DSP异构系统提出一种高效协同的矩阵乘加速方法，以提高矩阵乘的运算速度和最大化CPU+DSP异构系统的计算效率。技术方案是先初始化参数并对CPU+DSP异构系统信息配置，依据主处理器CPU和加速器DSP设计目标和计算性能的差异，将分配给计算结点待处理的数据划分给CPU和DSP协同处理，然后CPU和DSP并行进行数据传输与协同计算，得到个块矩阵C(i-1)(j-1)，最后将块矩阵C(i-1)(j-1)归并，组成M×N的结果矩阵C。采用本发明可使得CPU在负责数据传输和程序控制的同时积极与DSP协同完成矩阵乘计算，且数据传输与协同计算重叠，提高了CPU+DSP异构系统的矩阵乘运算速度和计算资源利用率。
【专利说明】面向CPU+DSP异构系统的矩阵乘加速方法
【技术领域】
[0001] 本发明涉及矩阵乘加速方法，尤指面向CPU+DSP异构计算系统的矩阵乘加速方 法。
【背景技术】
[0002] 异构计算系统是由主处理器和加速器两种不同体系结构处理器搭建的计算机系 统。目前，常见的异构计算系统有CPU+GPU和CPU+MIC组成异构计算系统，随着通用DSP计 算性能的不断攀升和通用DSP的广泛应用，CPU+DSP必将成为异构计算系统重要的发展方 向。
[0003] 矩阵乘是数值计算中最常用的一类操作，很多应用中都包含矩阵乘的计算过 程，面向CPU+DSP异构系统设计高效的矩阵乘方法可以有效提高应用的计算速度和提升 CPU+DSP异构系统的计算效率，以达到CPU+DSP异构系统的设计目标。
[0004] 矩阵乘是将被乘矩阵A的一行和乘数矩阵B的一列相乘得到结果矩阵C中的一个 元素。面向异构系统的矩阵乘法通常需要将矩阵乘计算过程合理分布在主处理器和加速器 之间协同完成计算过程，以提高矩阵乘的运算速度和最大化异构系统的计算效率。
[0005] 由于主处理器（CPU)和加速器（DSP或MIC或GPU)设计目标和计算性能的不同， 传统的面向异构系统的矩阵乘法通常由主处理器负责数据传输和程序控制，加速器负责矩 阵乘计算。该方法虽然也能够提高矩阵乘的运算速度，但是主处理器在加速器忙于计算的 过程中却处于等待计算结果的空闲状态，忽略了主处理器也拥有较高的计算性能，严重影 响了整个异构系统的计算效率。为了充分利用加速器和主处理器的计算潜能，最大化异构 计算系统的计算效率，E. Phillips等人率先在CPU+GPU异构系统上设计实现了协同的矩阵 乘，该方法中首先完成矩阵数据在CPU和GPU之间的划分，然后由主处理器CPU负责数据的 传输和程序控制，同时参与矩阵乘计算。这种高效协同的矩阵乘计算方法有效提高了矩阵 乘的运算速度和提升了 CPU+GPU异构系统的计算效率，达到了异构系统的设计目标。但是， 这种数据划分和传输方法以及控制策略仅适用于NVIDIA的统一架构GPU平台。
[0006] 综上所述，传统的异构系统矩阵乘方法无法满足CPU+DSP异构计算系统的设计目 标。因此，面向如图1所示的CPU+DSP异构系统，提供高效协同的矩阵乘方法加速方法是本 领域技术人员亟需解决的技术难题。
【发明内容】
[0007] 本发明要解决的技术问题在于：面向CPU+DSP异构系统提出一种高效协同的矩阵 乘加速方法，该方法中CPU在负责数据传输和程序控制的同时还积极与DSP协同完成矩阵 乘计算，且数据传输与协同计算重叠，以提高矩阵乘的运算速度和最大化CPU+DSP异构系 统的计算效率。
[0008] 为了解决上述技术问题，本发明的具体技术方案为：
[0009] 第一步、初始化参数并对CPU+DSP异构系统信息配置。具体步骤如下：
[0010] I. 1定义矩阵A的维度为MXK，矩阵B的维度为KXN，则A与B相乘的结果矩阵C 的维度为MXN，M，K，N均为正整数；
[0011] 1. 2查询CPU+DSP异构计算系统体系结构文档获取计算节点配置，即一个计算节 点由C颗CPU和D颗DSP组成；
[0012] 1. 3查询CPU+DSP异构计算系统体系结构文档获取计算结点的拓扑结构mXn，即， CPU+DSP异构系统由mXn个计算节点组成，每行有n个计算节点，每列有m个计算节点；
[0013] L 4依据计算结点中DSP的个数，将DSP分别标识为0, 1，d，（D-I)，0彡d < D ;
[0014] 1. 5依据CPU+DSP异构计算系统提供的基础软件工具集中的相关函数（如init) 完成DSP初始化；
[0015] 1. 6查询CPU+DSP异构计算系统体系结构文档，获取各处理器体系结构信息，即每 个计算结点中用于数据传输和流程控制的CPU核个数pd、用于计算的CPU核个数prt、每个 用于计算的CPU核拥有的浮点乘累加功能部件的数目m及主频f，每个DSP单元拥有的浮点 乘累加功能部件的数目m'及主频f ;
【权利要求】
1. 一种面向CPU+DSP异构系统的矩阵乘加速方法，其特征在于包括以下步骤： 第一步、初始化参数并对CPU+DSP异构系统信息配置，具体步骤如下： I. 1定义矩阵A的维度为MXK，矩阵B的维度为KXN，则A与B相乘的结果矩阵C的维 度为MXN，M，K，N均为正整数； 1. 2查询CPU+DSP异构计算系统体系结构文档获取计算节点配置，即一个计算节点由C 颗CPU和D颗DSP组成； 1. 3查询CPU+DSP异构计算系统体系结构文档获取计算结点的拓扑结构mXn，即， CPU+DSP异构系统由mXη个计算节点组成，每行有η个计算节点，每列有m个计算节点； 1. 4依据计算结点中DSP的个数，将DSP分别标识为0, 1，…d，…（D-I)，0彡d<D; 1. 5依据CPU+DSP异构计算系统提供的基础软件工具集中的相关函数完成DSP初始 化； 1. 6查询CPU+DSP异构计算系统体系结构文档，获取各处理器体系结构信息，即每个计 算结点中用于数据传输和流程控制的CPU核个数pd、用于计算的CPU核个数prt、每个用于 计算的CPU核拥有的浮点乘累加功能部件的数目m及主频f，每个DSP单元拥有的浮点乘累 加功能部件的数目m'及主频f; 1. 7计算主处理器的理论计算峰值Rpeak =Pet; 1. 8计算加速器的理论计算峰值R'peak =D*m' *f' ； 1. 9确定矩阵数据划分因子η=R'peak/Rpeak ; 第二步、依据主处理器(PU和加速器DSP设计目标和计算性能的差异，将分配给计算结 点待处理的数据划分给CPU和DSP协同处理，具体方法如下： 2. 1将M*K的矩阵A划分成m*K的块矩阵Ai, 丨，「]表示上取整，具体分块方法 如下： 2.I. 1 令i= 1 ; 2. 1. 2选取矩阵A的第（i-1) *m+l行至第i*m行组成块矩阵Ai ; 2.I. 3i=i+1 ； Γμ- 一 2. 1.4若?<一，转2.I. 2,否则，转2.I. 5,进行尾部处理； Im 2. 1. 5选取矩阵A的第（i-l) *m+l行至第M行组成块矩阵Ai ; jul 2. 2将K*N的矩阵B划分成K*n的子块矩阵Bj, 1表示上取整,具体分块 η 方法如下： 2· 2· 1 令j= 1 ; 2. 2. 2选取矩阵B的第（j-1) *n列至第j*n列组成块矩阵Bj ; 2. 2. 3j=j+1 ； 2. 2. 4若./< .转2. 2. 2,否则，转2. 2. 5,进行尾部处理； 2. 2. 5选取矩阵B的第（j-1) *n列至第N列组成块矩阵Bj; 2. 3依据矩阵数据划分因子Il将块矩阵h划分为％和％两个子矩阵，具体方法如 下： 2. 3. 1 令j= 1 ; 2. 3. 2依据矩阵数据划分因子η将块矩阵Bj划分为1Bj和2B j两个子矩阵，子矩阵1Bj 的维度为K=Knn2Bj的维度为K*n2,并且Ii1,n2满足公式⑴：
2. 3. 3j=j+1 ； 2. 3.4若./s^ ,转2. 3. 2,否则,转第三步； 第三步、CPU和DSP并行进行数据传输与协同计算，具体方法如下： 3. 1 令i= 1，j= 1 ; 3. 2采用CPU+DSP异构系统数据传输函数将块矩阵数据Ai和2BJ专输至DSP端存储空 间，具体步骤如下： 3. 2. 1在DSP端申请大小为size=MXkXsizeof(aij)的存储空间，sizeof(aij)表示 矩阵A中元素a、的存储长度，单位为字节； 3. 2. 2采用CPU+DSP异构系统提供的数据传输函数将块矩阵Ai传输至DSP端存储空间， 将DSP端存储空间中存储块矩阵Ai的空间称为Da,; 3. 2. 3i=i+1 ； 3· 2· 4 在DSP端申请大小为size=k*n2Xsizeof(IDij)的存储空间，sizeof(IDij)表示 矩阵B中元素的存储长度为多少字节； 3. 2. 5采用CPU+DSP异构系统提供的数据传输函数将子矩阵2Bj传输至DSP端存储空 间，将DSP端存储空间中存储块矩阵％的空间称为； 3. 2. 6j=j+1 ； 3.3采用3.2. 1和3.2. 2所述方法传输块矩阵数据Ai传输至DAi: 3. 4采用3. 2. 4和3. 2. 5所述方法传输块子矩阵数据2BJ专输至； 3. 5CPU和DSP同时并行执行以下操作： 3. 5. 1用于数据传输和流程控制的CPU核上的主线程负责CPU和DSP之间的通信和交 互，同时，主线程创建两个子线程T。和Td分别控制CPU端和DSP端的矩阵计算； 3. 5. 2T。调用矩阵乘库函数完成X的矩阵计算，其计算结果为结果矩阵C的 土夬矩阵C(H) (η)的子矩阵1C(H) (η); 3. 5. 3Td调用面向DSP体系结构的矩阵库函数完成X2B(k的矩阵计算，其计算结 果为结果矩阵C的块矩阵_的子矩阵2C^imM; 3. 5. 4主线程将子矩阵2C^ _传回至CPU端存储空间； 3. 5.5 释放D 3. 5. 6由子矩阵1_组成块矩阵_第1至第Ii1列，由子矩阵2C^ _组成 块矩阵C(H)chl)第h+l至第N列，其中，1^+? =N; 3. 6j=j+1 ； 3. 7如果P转3. 4,否则，转3. 8,进行尾部计算； 3.8释放^\「· 3. 9i=i+1 ； 3. 10如果，转3. 3,否则，转第四步； 第四步、将个块矩阵_归并，组成MXN的结果矩阵C，具体方法如下： 4. 1 令i= 1，j= 1 ; 4.2由块矩阵(^_1)+1)的第1至第111行的第1列至第11列组成结果矩阵(：的第（1-1)* 111+1 行至i*m行的第（j-1) *n+l列至第j*n列； 4. 3j=j+1 ； 4. 4如果_/<f,转4. 2,否则，转4. 5,进行列尾部处理； Il 4. 5由块矩阵CiM_的第1至第m行的第1列至第N- (i-1) *n列组成结果矩阵C的 第（i-1) *m+l行至第i*m行的第（j-1) *n+l列至第N列； 4. 6i=i+1 ； 4. 7如果/< 一，转4. 2,否则,转4. 8,进行行尾部处理； 綱 4. 8由块矩阵_的第1行第M-(i-1) *m行的第1列至第N-(i-1) *n列组成结果 矩阵C的第（i-1) *m+l行至第M行的第（j-1) *n+l列至第N列； 第五步、结束。
【文档编号】G06F17/16GK104317768SQ201410544785
【公开日】2015年1月28日 申请日期:2014年10月15日 优先权日:2014年10月15日
【发明者】刘杰, 迟利华, 甘新标, 晏益慧, 徐涵, 胡庆丰, 蒋杰, 李胜国, 王庆林, 皇甫永硕, 崔显涛, 周陈 申请人:中国人民解放军国防科学技术大学