CN103440121A

CN103440121A - 一种面向向量处理器的三角矩阵乘法向量化方法

Info

Publication number: CN103440121A
Application number: CN2013103642797A
Authority: CN
Inventors: 刘仲; 窦强; 郭阳; 邓让钰; 田希; 龚国辉; 唐遇星; 王永文; 彭元喜; 陈海燕; 胡封林; 陈磊
Original assignee: National University of Defense Technology
Current assignee: National University of Defense Technology
Priority date: 2013-08-20
Filing date: 2013-08-20
Publication date: 2013-12-11
Anticipated expiration: 2033-08-20
Also published as: CN103440121B

Abstract

本发明公开了一种面向向量处理器的三角矩阵乘法向量化方法，其步骤为：步骤1：将被乘数三角矩阵T中的三角矩阵元素按行连续存储；步骤2：根据向量处理器的向量处理单元个数和向量处理单元的MAC部件个数对乘数矩阵B按列划分为若干个子矩阵B_i；步骤3：依次实现被乘数三角矩阵T与子矩阵B_i的乘法，计算结果存储在原子矩阵B_i的存储位置；步骤4：遍历完乘数矩阵的全部子矩阵B_i，判断是否还有未计算的子矩阵B_i，若有，更新i=i+1，转步骤3；若无，则执行步骤5；步骤5：三角矩阵乘法的计算完成。本发明具有原理简单、操作简便、能充分发挥向量处理器运算效率等优点。

Description

一种面向向量处理器的三角矩阵乘法向量化方法

技术领域

本发明主要涉及到三角矩阵乘法计算领域，特指一种面向向量处理器的三角矩阵乘法向量化方法。

背景技术

基本线性代数函数库(Basic Linear Algebra Subprograms，BLAS)是各种科学计算广泛引用的核心数学库，也是许多数学软件的核心组件。目前，BLAS广泛应用于物理学、电子工程学、生物学、经济学、计算科学等科学与工程计算。其中，三角矩阵与矩阵乘法（Triangularmatrix-matrix multiplication）和三角矩阵求解（Triangular matrix solve）两类函数库是最常被调用的函数库之一。由于三角矩阵的转置和逆仍然是三角矩阵，这两类函数库最终可归结于三角矩阵乘法计算：C=TB，其中T为上三角矩阵或下三角矩阵，B为普通矩阵，C为结果矩阵。因为三角矩阵乘法是计算密集和访存密集的运算，对处理器的运算能力、访存带宽及延迟的要求非常高，因此，针对特定体系结构的特点对三角矩阵乘法进行优化设计，对提高程序的运行速度、发挥处理器的运算能力具有很重要的意义。

另一方面，随着大型稠密线性方程组求解、雷达信号处理、高清视频和数字图像处理等计算密集型应用的高性能计算需求日益增长，计算机体系结构出现显著变化，出现许多新型体系结构。向量处理器体系结构是其中的一种新颖体系结构。如图1所示，向量处理器一般包括标量处理部件（Scalar Processing Unit，SPU）和向量处理部件（Vector Processing Unit，VPU），SPU负责标量任务计算和流控；VPU负责向量计算，包括若干向量处理单元（VectorProcessing Element，VPE），每个VPE包含MAC、ALU、BP等多个功能部件。SPU和VPU之间提供数据传输和交换数据的通道，向量处理器提供广播指令将SPU的一个标量寄存器的数据广播到VPU的一个向量寄存器。向量数据访问单元支持向量数据的Load/Store，提供大容量的专用向量存储器。向量处理器具有强大的计算能力，能够大幅度提高程序的计算性能，但同时对软件开发提出了新的挑战，如何针对向量处理器多处理单元、多功能部件等体系结构特点，充分开发各个层次的并行性，将现有应用高效地向量化是当前面临的主要困难。

对于普通的矩阵乘法C=A*B，其中A为m*n阶矩阵，B为n*k阶矩阵，C为m*k阶矩阵。结果矩阵C按照公式（1）计算：

C_{ij} = Σ_{r = 0}^{n - 1} A_{ir} * B_{rj} (0 \leq i \leq m - 1,0 \leq j \leq k - 1) - - - (1)

根据公式（1）计算结果矩阵C，在计算时间方面，每个C矩阵元素的计算需要n次乘法和n-1次加法，完成矩阵C的计算需要m*k*n次乘法和m*k*(n-1)次加法。在存储方面，若以矩阵元素为一个存储单位，则需要存储空间为（m*n+n*k+m*k）单位。与普通矩阵相比，三角矩阵有一半的矩阵元素为0，显然，简单的套用矩阵乘法实现三角矩阵乘法，在计算效率和存储效率方面都是低效的，需要针对三角矩阵计算和存储特点，以及所在处理器平台的体系结构特点制定更高效的计算方法。

目前，提高三角矩阵乘法计算性能的方法包括硬件方法和软件方法，如基于多FPGA的矩阵乘法并行计算系统（申请号：200710069953.3），一种基于龙芯3A的单精度矩阵乘法优化方法（申请号：201110191133.8），具有减少的带宽要求的矩阵乘法（申请号：200710097456.4）等。硬件方法速度快，但是成本高，不够灵活。软件的方法往往计算效率较低，且大多是矩阵与矩阵的乘法，没有专门考虑三角矩阵与矩阵的乘法。与普通矩阵相比，三角矩阵有一半的矩阵元素为0，因此在计算和存储上套用普通的矩阵乘法不能有效提升三角矩阵的计算性能。并且，在面向向量处理器平台上的计算，一方面由于三角矩阵不像普通矩阵是方阵，难以实现规整化的数据对齐，不能平衡各个处理单元的计算负载；另一方面，向量数据访问单元支持的向量数据Load/Store，一般只支持按行的向量数据读取，不支持按列的向量数据读取，不能支持如专利“具有减少的带宽要求的矩阵乘法”（申请号：200710097456.4）支持的按列访问矩阵元素。如何在向量处理器上高效实现三角矩阵乘法的向量化计算，才能充分发挥向量处理器的计算能力是目前需要解决的技术难点，在公开文献中也没有面向向量处理器的三角矩阵乘法向量化方法。

因此，需要设计一种面向向量处理器的高效的三角矩阵乘法向量化方法，以便最优的发挥向量处理器的运算效率。

发明内容

本发明要解决的技术问题就在于：针对现有技术存在的技术问题，本发明提供一种原理简单、操作简便、能充分发挥向量处理器运算效率的面向向量处理器的三角矩阵乘法向量化方法，从而充分利用向量处理器多处理单元、多功能部件的体系结构特点，开发向量处理器的指令、数据、任务等各个层次的并行性，实现了向量处理器的计算负载平衡。

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

一种面向向量处理器的三角矩阵乘法向量化方法，其步骤为：

步骤1：将被乘数三角矩阵T中的三角矩阵元素按行连续存储；设被乘数三角矩阵T为m*m阶矩阵，乘数矩阵B为m*k阶矩阵，计算三角矩阵T与乘数矩阵B的乘法，结果矩阵C=T*B，C为m*k阶矩阵；

步骤2：根据向量处理器的向量处理单元个数和向量处理单元的MAC部件个数对乘数矩阵B按列划分为若干个子矩阵B_i；

步骤3：依次实现被乘数三角矩阵T与子矩阵B_i的乘法，计算结果存储在原子矩阵B_i的存储位置；

步骤4：遍历完乘数矩阵的全部子矩阵B_i，判断是否还有未计算的子矩阵B_i，若有，更新i=i+1，转步骤3；若无，则执行步骤5；

步骤5：三角矩阵乘法的计算完成。

作为本发明的进一步改进：

所述步骤（2）中，设向量处理器的向量处理单元个数为p，向量处理单元的MAC部件个数为q；对乘数矩阵B按列划分为子矩阵，子矩阵的行数与B矩阵一致，均为m，子矩阵的列数固定为p*q，若k不是(p*q)的整数倍，则最后一个子矩阵的列数为k除以(p*q)的余数；记子矩阵的个数为s，子矩阵依次记为B₀,B₁,…,B_s-1；采用分块算法实现结果矩阵C的计算：(C₀,C₁,…,C_s-1)=T*(B₀,B₁,…,B_s-1)，令i=0。

所述步骤（3）中，具体步骤为：

步骤3.1：设定初始参数，若三角矩阵T为上三角矩阵，则P指向三角矩阵T的第一行，Q和S分别指向子矩阵B_i的第一行；若三角矩阵T为下三角矩阵，则P指向三角矩阵T的最后一行，Q和S分别指向子矩阵B_i的最后一行；

步骤3.2：向量处理器的SPU用标量存取指令读取三角矩阵T的第P行中的连续q个元素到q个标量寄存器，并用向量处理器提供的广播指令分别广播到VPU的q个向量寄存器；

步骤3.3：向量处理器的VPU用向量存取指令读取子矩阵B_i的第S行的q*p个元素到与上述步骤3.2不同的q个向量寄存器；

步骤3.4：向量处理器的VPU对上述步骤3.2的q个向量寄存器和步骤3.3的q个向量寄存器分别执行乘法，乘法结果分别累加到与上述步骤3.2和3.3不同的q个向量寄存器；

步骤3.5：判断子矩阵B_i中是否还有另一行数据，若有，则子矩阵B_i中更新S为S的下一行（三角矩阵T为上三角矩阵情况）或上一行（三角矩阵T为下三角矩阵情况），转步骤3.2；

步骤3.6：将上述步骤3.4计算结果的q个向量寄存器保存到子矩阵B_i的第Q行；

步骤3.7：判断三角矩阵T中是否还有另一行数据，若有，则三角矩阵T中更新P为P的下一行（三角矩阵T为上三角矩阵情况）或上一行（三角矩阵T为下三角矩阵情况）；子矩阵B_i中更新Q为Q的下一行（三角矩阵T为上三角矩阵情况）或上一行（三角矩阵T为下三角矩阵情况），子矩阵B_i中更新S为Q行，转步骤3.2；

步骤3.8；被乘数三角矩阵T与子矩阵B_i的乘法计算完成。

与现有技术相比，本发明的优点在于：

（1）本发明充分利用三角矩阵中有一半元素为0，0元素不需要与矩阵B中的对应行的数据相乘，这样的方法能够大幅减少计算量。

（2）本发明中被乘数矩阵T中的三角矩阵元素按行连续存储，不参与计算的0元素不需要存储；实现了原位计算，结果矩阵是保存在乘数矩阵B中，节省了存储空间。

（3）本发明中，每次向量计算都能够平衡各个处理单元的计算负载，充分发挥向量处理器的各个处理单元的计算能力。

（4）本发明中，避免了通常向量处理器不支持的矩阵列向量数据的访问和向量处理单元之间的浮点归约求和（浮点归约求和开销很大）。本发明实现了高效的三角矩阵乘法的向量化计算，从而大幅度提高了三角矩阵乘法的计算性能。

附图说明

图1是向量处理器的一般结构原理示意图。

图2是本发明方法的流程示意图。

图3是本发明方法的被乘数三角矩阵T与子矩阵B_i的乘法计算的流程示意图。

图4是本发明中被乘数上三角矩阵T与子矩阵B_i乘法在具体实施例中的运算示意图。

具体实施方式

以下将结合说明书附图和具体实施例对本发明做进一步详细说明。

如图2所示，本发明面向向量处理器的三角矩阵乘法向量化方法，具体步骤为：

步骤2：根据向量处理器的向量处理单元个数和向量处理单元的MAC部件个数对乘数矩阵B按列划分为若干个子矩阵B_i；设向量处理器的向量处理单元个数为p，向量处理单元的MAC部件个数为q。对乘数矩阵B按列划分为子矩阵，子矩阵的行数与B矩阵一致，均为m，子矩阵的列数固定为p*q，若k不是(p*q)的整数倍，则最后一个子矩阵的列数为k除以(p*q)的余数。记子矩阵的个数为s，子矩阵依次记为B₀,B₁,…,B_s-1。采用分块算法实现结果矩阵C的计算：(C₀,C₁,…,C_s-1)=T*(B₀,B₁,…,B_s-1)，令i=0；

步骤3：依次实现被乘数三角矩阵T与子矩阵B_i的乘法，计算结果存储在原子矩阵B_i的存储位置；计算被乘数三角矩阵T与子矩阵B_i的乘法：C_i=T*B_i，计算结果存储在原子矩阵B_i的存储位置；

步骤4：遍历完乘数矩阵的全部子矩阵B_i；判断是否还有未计算的子矩阵B_i，若有，更新i=i+1，转步骤3；若无，则执行步骤5；

步骤5：三角矩阵乘法的计算完成。

如图3所示，本实施例的上述步骤3中，被乘数三角矩阵T与子矩阵B_i的乘法计算：C_i=T*B_i，其详细步骤为：

步骤3.1：设定初始参数，若三角矩阵T为上三角矩阵，则P指向三角矩阵T的第一行，Q和S分别指向子矩阵B_i的第一行。若三角矩阵T为下三角矩阵，则P指向三角矩阵T的最后一行，Q和S分别指向子矩阵B_i的最后一行；

步骤3.7：判断三角矩阵T中是否还有另一行数据，若有，则三角矩阵T中更新P为P的下一行（三角矩阵T为上三角矩阵情况）或上一行（三角矩阵T为下三角矩阵情况）。子矩阵B_i中更新Q为Q的下一行（三角矩阵T为上三角矩阵情况）或上一行（三角矩阵T为下三角矩阵情况），子矩阵B_i中更新S为Q行，转步骤3.2；

步骤3.8；被乘数三角矩阵T与子矩阵B_i的乘法计算完成。

在上述方法中，被乘数矩阵T中的三角矩阵元素按行连续存储，不参与计算的0元素不需要存储；结果矩阵是保存在乘数矩阵B中，实现了原位计算。若设T为m*m阶矩阵，B为m*k阶矩阵，则存储空间为（m*(m+1）/2+m*k）元素，相比普通矩阵乘法减少了（m*(m-1）/2+m*k）元素。其原理是充分利用了三角矩阵中有一半元素为0，0元素不需要与矩阵B中的对应行的数据相乘。这样的方法既减少了计算量，又节省了存储空间，并且充分发挥了向量处理器的各个处理单元的计算能力，平衡了各个处理单元的计算负载，实现了高效的三角矩阵乘法的向量化计算，从而大幅度提高了三角矩阵乘法的计算性能。

如图4所示，为本发明中被乘数上三角矩阵T与子矩阵B_i乘法在具体实施例的运算示意图。

在本实例中，以一个上三角矩阵T与子矩阵B_i的乘法计算实例阐述本发明方法的实施过程。

如图4(a)所示，设向量处理器的VPE数目为4，VPE的MAC部件个数为1，被乘数上三角矩阵T的行数与列数等于4，子矩阵B_i的行数与列数等于4。实现被乘数上三角矩阵T与子矩阵B_i乘法的方法和步骤是：

（1）设定初始参数，因为实施例中三角矩阵T为上三角矩阵，设定P指向三角矩阵T的第一行，Q和S分别指向子矩阵B_i的第一行；即P=0,Q=0,S=0；

（2）如图4（b），向量处理器的SPU用标量存取指令读取三角矩阵T的第0行中的1个元素t00元素到标量寄存器，通过向量处理器的广播指令广播到1个向量寄存器(t00,t00,t00,t00)。向量处理器的VPU用向量存取指令读取子矩阵B_i的第0行的4个元素到1个向量寄存器(b00,b01,b02,b03)。向量处理器的VPU对上述2个向量寄存器执行乘法，乘法结果累加到1个向量寄存器。更新S=S+1，依次下去，当S遍历完0，1，2，3后，得到结果矩阵的第0行，并保存到子矩阵B_i的第0行B0。

（3）如图4（c），更新参数P、Q、S。即P=1,Q=1,S=1。向量处理器的SPU用标量存取指令读取三角矩阵T的第1行中的1个元素t11元素到标量寄存器，通过向量处理器的广播指令广播到1个向量寄存器(t11,t11,t11,t11)。向量处理器的VPU用向量存取指令读取子矩阵B_i的第1行的4个元素到1个向量寄存器(b10,b11,b12,b13)。向量处理器的VPU对上述2个向量寄存器执行乘法，乘法结果累加到1个向量寄存器。更新S=S+1，依次下去，当S遍历完1，2，3后，得到结果矩阵的第1行，并保存到子矩阵B_i的第1行B1。

（4）如图4（d），更新参数P、Q、S。即P=2,Q=2,S=2。向量处理器的SPU用标量存取指令读取三角矩阵T的第2行中的1个元素t22元素到标量寄存器，通过向量处理器的广播指令广播到1个向量寄存器(t22,t22,t22,t22)。向量处理器的VPU用向量存取指令读取子矩阵B_i的第2行的4个元素到1个向量寄存器(b20,b21,b22,b23)。向量处理器的VPU对上述2个向量寄存器执行乘法，乘法结果累加到1个向量寄存器。更新S=S+1，依次下去，当S遍历完2，3后，得到结果矩阵的第2行，并保存到子矩阵B_i的第2行B2。

（5）如图4（e），更新参数P、Q、S。即P=3,Q=2,S=3。向量处理器的SPU用标量存取指令读取三角矩阵T的第3行中的1个元素t33元素到标量寄存器，通过向量处理器的广播指令广播到1个向量寄存器(t33,t33,t33,t33)。向量处理器的VPU用向量存取指令读取子矩阵B_i的第3行的4个元素到1个向量寄存器(b30,b21,b32,b33)。向量处理器的VPU对上述2个向量寄存器执行乘法，乘法结果累加到1个向量寄存器。更新S=S+1，依次下去，当S遍历完3后，得到结果矩阵的第3行，并保存到子矩阵B_i的第3行B3。

（6）P遍历完三角矩阵的所有行0，1，2，3后，三角矩阵T与子矩阵B_i的乘法计算完成。

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

Claims

1.一种面向向量处理器的三角矩阵乘法向量化方法，其特征在于，步骤为：

步骤5：三角矩阵乘法的计算完成。

2.根据权利要求1所述的面向向量处理器的三角矩阵乘法向量化方法，其特征在于，所述步骤（2）中，设向量处理器的向量处理单元个数为p，向量处理单元的MAC部件个数为q；对乘数矩阵B按列划分为子矩阵，子矩阵的行数与B矩阵一致，均为m，子矩阵的列数固定为p*q，若k不是(p*q)的整数倍，则最后一个子矩阵的列数为k除以(p*q)的余数；记子矩阵的个数为s，子矩阵依次记为B₀,B₁,…,B_s-1；采用分块算法实现结果矩阵C的计算：(C₀,C₁,…,C_s-1)=T*(B₀,B₁,…,B_s-1)，令i=0。

3.根据权利要求1所述的面向向量处理器的三角矩阵乘法向量化方法，其特征在于，所述步骤（3）中，具体步骤为：

步骤3.5：判断子矩阵B_i中是否还有另一行数据，若有，则子矩阵B_i中更新S为S的下一行，即三角矩阵T为上三角矩阵情况；或上一行，即三角矩阵T为下三角矩阵情况；转步骤3.2；

步骤3.7：判断三角矩阵T中是否还有另一行数据，若有，则三角矩阵T中更新P为P的下一行，即三角矩阵T为上三角矩阵情况；或上一行，即三角矩阵T为下三角矩阵情况；子矩阵B_i中更新Q为Q的下一行，即三角矩阵T为上三角矩阵情况；或上一行，即三角矩阵T为下三角矩阵情况，子矩阵B_i中更新S为Q行，转步骤3.2；

步骤3.8；被乘数三角矩阵T与子矩阵B_i的乘法计算完成。