[发明专利]面向异构融合体系结构的矩阵乘加速方法

权利要求书

1.一种面向异构融合体系结构的矩阵乘加速方法，其特征在于包括以下步骤：

第一步、设计面向异构融合体系结构的分块矩阵乘版本，具体步骤如下：

1.1异构融合体系结构配置与初始化，具体方法是：

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

1.2若异构融合体系结构仅由CPU组成，对CPU进行初始化，方法如下：

1.2.1.查询体系结构手册获取一个CPU中计算核的拓扑结构，如果多核CPU中的拓扑结构为mg×ng，即一个CPU中有mg×ng个计算核，物理分布为mg行×ng列，将计算核依次编号为(0,0),(0,1),…(0,ng-1)，(1,0),(1,1),…(1,ng-1)，……，(i,0),(i,1),…,(i,j),(i,ng-1)，……(mg-1,0),(mg-1,1),…(mg-1,ng-1)，0≤i≤mg-1，0≤j≤ng-1；

1.2.2.查询体系结构手册获取每个CPU计算核拥有的浮点向量乘累加功能部件FMA的数目m_e；

1.2.3.利用操作系统提供的初始化函数完成多核CPU初始化；

1.3依据CPU的拓扑结构对A和B进行矩阵划分，将矩阵A划分为mg×ng个A块矩阵，mg×ng个A块矩阵的排布方式与mg×ng个计算核的排布方式相同，第i行第j列个A块矩阵用A_ij表示，0≤i≤mg-1,0≤j≤ng-1，每个块矩阵的维度为m×k，其中表示上取整；将K*N的矩阵B划分为ng×mg个B块矩阵，ng×mg个B块矩阵的排布方式与mg×ng个计算核的排布方式相同，第r行第s列B块矩阵用B_rs表示，0≤r≤mg-1,0≤s≤ng-1，每个块矩阵的维度为k×n，

1.4初始化结果矩阵C＝0，即将C中各元素均赋值为0，令C(mo,no)＝0，C(mo,no) 表示结果矩阵第mo行第no列的元素，0≤mo≤M-1,0≤no≤N-1；

1.5矩阵乘加速，得到mg×ng×m_e个C块矩阵：mg×ng×m_e个FMA以CPU计算核共享存储空间为共享存储空间，以标量数据缓冲区为数据缓冲区，并行执行分块矩阵乘法操作，使得矩阵乘加速，每个FMA独立完成分配给自己的块矩阵乘，编号为(i',j')的FMA完成块矩阵乘

1.6结果归并，即依据数据分布原理，以CPU计算核共享存储空间为共享存储空间，将共享存储空间的结果归并组成结果矩阵C，按1.2～1.6采用MPI编程方式实现的矩阵乘即为CPU矩阵乘加速版本v_cpu；

1.7若异构融合体系结构由CPU和GPU组成，采用CUDA编程方式实现面向GPU的矩阵乘加速版本v_gpu,具体方法如下：

1.7.1查询体系结构手册获取每个GPU中计算核的拓扑结构，即GPU中的线程处理器TPU拓扑结构mg×ng，即GPU中有mg×ng个TPU，物理分布为mg行×ng列，依次编号为(0,0),(0,1),…(0,ng-1)，

(1,0),(1,1),…(1,ng-1)，……，(i,0),(i,1),…,(i,j),…,(i,ng-1),……(mg-1,0),(mg-1,1),…(mg-1,ng-1)；

1.7.2查询体系结构手册获取每个TPU拥有的向量浮点乘累加功能部件FMA的数目m_e；

1.7.3利用操作系统提供的函数完成初始化；

1.7.4矩阵划分，基于GPU的拓扑结构采用1.3步所述的矩阵划分方法对A和B进行分块，得到mg×ng个A块矩阵和ng×mg个B块矩阵；

1.7.5 mg×ng×m_e个FMA以GPU的Shared Memory为共享存储空间，以GPU的ConstantMemory为数据缓冲区，采用1.5步所述的矩阵乘加速方法并行执行矩阵乘法，每个FMA独立完成分配给自己的块矩阵乘，编号为(i',j')的FMA完成块矩阵乘

1.7.6以GPU的Shared Memory为共享存储空间，采用1.6步所述的结果归并方法完成各FMA计算的C块矩阵结果归并，得到C；

1.8若异构融合体系结构由CPU与MIC组成，采用Offload编程方式实现的MIC矩阵乘加速版本v_mic,具体方法如下：

1.8.1查询体系结构手册获取每个MIC中计算核，即如MIC中的向量处理器VPU拓扑结构mg×ng，即一个MIC中有mg×ng个VPU，物理分布为mg行×ng列，依次编号为(0,0),(0,1),…(0,ng-1)，(1,0),(1,1),…(1,ng-1)，……，(i,0),(i,1),…,(i,j),…,(i,ng-1),……(mg-1,0),(mg-1,1),…(mg-1,ng-1)；

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

1.8.3利用操作体系提供的函数完成初始化；

1.8.4矩阵划分，基于VPU的拓扑结构采用1.3步所述的矩阵划分方法对A和B进行分块，得到mg×ng个A块矩阵和ng×mg个B块矩阵；

1.8.5 mg×ng×m_e个FMA以MIC的内存空间即内存为共享存储空间，以共享存储空间便携式Cache为数据缓冲区，采用1.5步所述的矩阵乘加速方法并行执行矩阵乘法，每个FMA独立完成分配给自己的块矩阵乘，编号为(i',j')的FMA完成块矩阵乘

1.8.6以MIC的内存空间为共享存储空间，采用1.6步所述的结果归并方法完成各FMA计算的C块矩阵结果归并，得到C；

1.9若异构融合体系结构由CPU与Matrix2000组成，若Matrix2000采用SCIF编程模式，则采用SCIF编程模式实现矩阵乘加速版本v_scif,具体方法如下：

1.9.1查询体系结构手册获取每个Matrix2000中计算核，即Matrix2000中的向量处理单元VPE拓扑结构mg×ng，即Matrix2000中有mg×ng个VPE，物理分布为mg行×ng列，依次编号为(0,0),(0,1),…(0,ng-1)，(1,0),(1,1),…(1,ng-1)，……，(i,0),(i,1),…,(i,j),…,(i,ng-1),……(mg-1,0),(mg-1,1),…(mg-1,ng-1)；

1.9.2查询体系结构手册获取每个VPE拥有的向量浮点乘累加功能部件FMA(FusedMultiply Add)的数目m_e；

1.9.3利用操作系统提供的函数完成初始化；

1.9.4矩阵划分，基于VPE的拓扑结构采用1.3步所述的矩阵划分方法对A和B进行分块，得到mg×ng个A块矩阵和ng×mg个B块矩阵；

1.9.5 mg×ng×m_e个FMA以Matrix2000的全局空间Global Cache为共享存储空间，以共享存储空间Array Memory为数据缓冲区，采用1.5步所述的矩阵乘加速方法并行执行矩阵乘法；每个FMA独立完成分配给自己的块矩阵乘，编号为(i',j')的FMA完成块矩阵乘

1.9.6以Matrix2000中的共享存储空间Array Memory为共享存储空间，采用1.6步所述的结果归并方法完成各FMA计算的C块矩阵结果归并，得到C；

1.10若异构融合体系结构由CPU与Matrix2000组成，若Matrix2000采用COI编程模式，则采用COI编程模式实现矩阵乘加速版本v_coi,具体方法如下：

1.10.1查询体系结构手册获取系统中计算核，若Matrix2000中的向量处理单元VPE拓扑结构为mg×ng，即Matrix2000中有mg×ng个VPE，物理分布为mg行×ng列，依次编号为(0,0),(0,1),…(0,ng-1)，(1,0),(1,1),…(1,ng-1)，……，(i,0),(i,1),…,(i,j),…,(i,ng-1),……(mg-1,0),(mg-1,1),…(mg-1,ng-1)；

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

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

1.10.4矩阵划分，基于VPE的拓扑结构采用1.3步所述的矩阵划分方法对A和B 进行分块，得到mg×ng个A块矩阵和ng×mg个B块矩阵；

1.10.5 mg×ng×m_e个FMA以Matrix2000的全局空间Global Cache为共享存储空间，以共享存储空间Array Memory为数据缓冲区，采用1.5步所述的矩阵乘加速方法并行执行矩阵乘法；每个FMA独立完成分配给自己的块矩阵乘，编号为(i',j')的FMA完成块矩阵乘

1.10.6以Matrix2000中的共享存储空间Array Memory为共享存储空间，采用1.6步所述的结果归并方法完成各FMA计算的C块矩阵结果归并，得到C；

1.11若异构融合体系结构由CPU与Matrix2000组成，若Matrix2000采用OpenMP taget编程模式，则采用OpenMP taget编程模式实现OpenMP taget矩阵乘加速版本v_target,具体方法如下：

1.11.1查询体系结构手册获取系统中计算核，即如Matrix2000中的向量处理单元VPE拓扑结构mg×ng，即系统中有mg×ng个VPE，物理分布为mg行×ng列，依次编号为(0,0),(0,1),…(0,ng-1)，(1,0),(1,1),…(1,ng-1)，……，(i,0),(i,1),…,(i,j),…,(i,ng-1),……(mg-1,0),(mg-1,1),…(mg-1,ng-1)；

1.11.2查询体系结构手册获取每个VPE拥有的向量浮点乘累加功能部件FMA(FusedMultiply Add)的数目m_e；

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

1.11.4矩阵划分，基于VPE的拓扑结构采用1.3步所述的矩阵划分方法对A和B进行分块，得到mg×ng个A块矩阵和ng×mg个B块矩阵；

1.11.5基于Matrix2000的全局空间Global Cache和共享存储空间Array Memory，mg×ng×m_e个FMA以Matrix2000的全局空间Global Cache为共享存储空间，以共享存储空间Array Memory为数据缓冲区，采用1.5步所述的矩阵乘加速方法并行执行矩阵乘法；每个FMA独立完成分配给自己的块矩阵乘，编号为(i',j')的FMA完成块矩阵乘

1.11.6以Matrix2000中的共享存储空间Array Memory为共享存储空间，采用1.6步所述的结果归并方法完成各FMA完成的C块矩阵结果归并，得到C；

第二步：集成异构融合多版本矩阵乘，方法如下：

2.1分别编译v_cpu、v_gpu、v_mic、v_target、v_coi和v_scif对应的源代码，生成v_cpu、v_gpu、v_mic、v_target、v_coi和v_scif的可执行文件；

2.2利用tar命令将v_cpu、v_gpu、v_mic、v_target、v_coi和v_scif的可执行文件打包生成异构融合版本的库文件HU-xgemm；

第三步、适配异构融合体系结构中的加速器，方法如下：

3.1利用操作系统命令lspci|grep processor查询异构融合体系结构中加速器类型，令加速器类型变量arc为加速器类型；

3.2若arc＝Matrix2000，转3.2.1：否则，转3.3；

3.2.1.查询Matrix2000编程技术手册，确认支持的编程模型，赋值给编程模型变量prolan；

3.2.2.若prolan＝OpenMP taget，调用HU-xgemm中的v_cpu和v_target分别完成主处理器CPU端和加速器Matrix2000端的矩阵乘计算；否则，转3.2.3；

3.2.3.若prolan＝COI，调用HU-xgemm中的v_cpu和v_coi分别完成主处理器CPU端和加速器Matrix2000端的矩阵乘计算；否则，转3.2.4；

3.2.4.若prolan＝SCIF，调用HU-xgemm中的v_cpu和v_scif分别完成主处理器CPU端和加速器Matrix2000端的矩阵乘计算；否则，转3.3；

3.3若arc＝MIC，调用HU-xgemm中的v_cpu和v_mic分别完成主处理器CPU端和加速器MIC端的矩阵乘计算；否则，转3.4；

3.4若arc＝GPU，调用HU-xgemm中的v_cpu和v_gpu分别完成主处理器CPU端和加速器GPU端的矩阵乘计算；否则，转3.5；

3.5若系统中无专用加速器，即只有CPU，调用HU-xgemm中的v_cpu完成矩阵乘计算。

2.如权利要求1所述的面向异构融合体系结构的矩阵乘加速方法，其特征在于1.2.3步所述操作系统提供的初始化函数指init函数。