CN102446160A - 面向双精度simd部件的矩阵乘实现方法 - Google Patents

Abstract

本发明公开了一种面向双精度SIMD部件的矩阵乘实现方法，目的是提高矩阵乘在SIMD部件上的计算速度。技术方案是先增加矩阵B和矩阵C的列数；然后对矩阵A、B进行分块；A的每个分块和对应的B分块使用SIMD部件进行相乘，并将结果加到矩阵C中相应位置的结果和上；采用本发明避免了对矩阵数据进行重排序，提高了矩阵乘在SIMD部件上的计算速度。

Description

面向双精度SIMD部件的矩阵乘实现方法

技术领域

本发明涉及通用微处理器中的SIMD(单指令多数据)部件，尤其指面向双精度SIMD部件的矩阵乘实现方法。

背景技术

通用微处理器芯片的集成度越来越大，在处理器中设计实现支持双精度浮点计算的SIMD部件来支持大规模科学和工程计算是一个重要的发展趋势。目前，商用微处理器上已经集成了SIMD部件，如Intel的MMX/SSE/AVX以及AMD的3D Now！技术等，都是面向SIMD部件的SIMD指令集。SIMD部件利用SIMD指令对向量进行操作，一个向量由多个浮点数据组成，从而实现了同时对多个浮点数据进行操作，加速了计算过程。

矩阵乘操作是数值计算中最常用的一类操作，很多应用中都包含矩阵乘的计算过程，利用SIMD部件加速矩阵乘计算过程可以有效提高应用的计算速度。实现面向SIMD部件的高效矩阵乘方法是关系到发挥SIMD部件加速能力的关键。否则，难以达到SIMD部件的加速计算的设计目标。

矩阵乘法是将被乘数矩阵A的一行和乘数矩阵B的一列相乘得到结果矩阵C中的一个元素。由于访问A和B的顺序不同，要对矩阵A或B进行重排序，传统的做法是将其中的一个矩阵进行转置。中国专利200380107095.7提出了一种使用SIMD寄存器的小矩阵乘法，方法中首先对矩阵数据进行重排序及寄存器加载，随后，将被乘数矩阵A的对角线乘以乘数矩阵B的列，把结果被加到结果矩阵C的列的结果和上。但是，这种方法只能处理维度较小的两个矩阵相乘。2001年Douglas和Jonathan提出了一种面向Intel SIMD部件的矩阵乘实现方法，但是，这种方法只适用于Intel的SIMD部件，而且方法中首先要对输入矩阵B进行转置操作。美国专利US007873812B1提出了一种面向SIMD部件的矩阵乘实现方法，但是它对矩阵的列数有特殊要求，只能处理输入矩阵A的列数可以被W(SIMD部件的向量长度)整除的情况，而且需要先对输入矩阵A进行转置操作，并要使用选择部件来选择向量寄存器中的元素。

综上所述，这些方法中都要求对矩阵数据进行有效的重排序，重排序会导致较大的计算开销，影响了矩阵乘在SIMD部件上的计算速度。面向SIMD部件，实现无需重排序的矩阵相乘方法仍是本领域技术人员迫切希望解决的技术问题。

发明内容

本发明要解决的技术问题是提出一种面向包含广播指令ld1toW的双精度SIMD部件的矩阵乘实现方法，避免对矩阵数据进行重排序，提高矩阵乘在SIMD部件上的计算速度。广播指令ld1toW是指将1个双精度浮点数据从存储器装载到向量寄存器的W个位置中。

本发明的技术方案为：对矩阵A和B进行分块，使用SIMD对A、B的子矩阵相乘，将子矩阵结果累加到结果矩阵C。

具体的技术方案为：

设A和B为输入矩阵，且均为双精度浮点矩阵，结果矩阵为C(C＝A×B)，A矩阵的大小是M*K，B矩阵为K*N，C大小为M*N。双精度SIMD部件的向量长度为W，即一个向量包含W个双精度浮点数。N、W、K全是整数。

第一步，增加矩阵B的列数，对增加的列使用0进行数据填充；

SIMD部件每次按行读取矩阵B的W个数据，如果N不是W的整数倍，在对B中每行数据的最后一次读取时，不能得到和矩阵A相乘的正确数据，这样就会得到错误结果。

所以，当N不是W的整数倍时，增加B的列数，将B增加W-N％W列，使B的列数为N+W-N％W，％表示模运算，增加的列使用0进行数据填充；当N是W的整数倍时，B的列数不变。

第二步，增加矩阵C的列数并将矩阵C的内容初始化为全0；

矩阵乘结果使用向量存指令存储到矩阵C中，矩阵C的列数必须和矩阵B相同，因此，需要增加矩阵C的列数，使C的列数为N+W-N％W。

矩阵C需要存储计算的中间结果，并对中间结果进行累加，所以需要将矩阵C的初始值初始化为0。

第三步，根据SIMD部件的向量寄存器数目VN对矩阵B进行分块，将K*N的矩阵B划分成k*n的子块B_j

其中n必须是W的整数倍；VN为正整数。

当矩阵B的K*N较大时，SIMD部件不可能将B中的所有数据取至SIMD部件的寄存器中，为了提高计算效率，需要对矩阵B进行分块。这样可以使每次子矩阵相乘的过程中，矩阵B的数据在SIMD部件的寄存器中被重复利用，提高SIMD部件的计算效率。

子矩阵B_j的大小即n和k值须满足：

2+n*(k+1)/W＜VN且n％W＝0和k％W＝0，

2+n*(k+1)/W为每次子矩阵相乘所需的最少的向量寄存器数目(使用1个向量寄存器存储矩阵A的一个数据，使用n*k/W个向量寄存器存储B矩阵的数据，使用1个向量寄存器存储向量乘结果，使用n/W个向量寄存器存储每行的最终计算结果)。

第四步，将M*K的矩阵A划分成M*k的子块A_i

表示下取整。

第五步，子矩阵A_i和B_j在SIMD部件中相乘，并将结果累加到结果矩阵C中；

5.1令i＝1，j＝1，u＝1，v＝1；

5.2将n/W个结果向量寄存器V_s的内容初始化为0，1≤s≤n/W；

5.3使用广播指令ld1toW将A_i中的一个元素a_uv取至向量寄存器V₀；

5.4令P＝1；

5.5如果u等于1，使用向量访存指令将B_j中的第v行元素中从(P-1)*W+1至P*W的元素取至向量寄存器V_Z中，1+n/W≤z≤n*(k+1)/W，执行第5.6步；如果u不等于1，则数据已经存放在向量寄存器中，执行第5.6步；

5.6V₀和V_Z进行向量乘法操作，将结果存储在向量寄存器V_t中t＝1+n*(k+1)/W；

5.7V_t和结果向量寄存器V_s进行向量加操作，将结果存放在V_s中；

5.8如果P＜n/W，P＝P+1，跳转至5.5；否则，执行5.9步；

5.9将n/W个结果向量寄存器V_s中的数据和C的第u行中第(u-1)*n+1+(i-1)*n列至第u*n+(i-1)*n列的n个数据进行累加，并将结果写到C中；

5.10如果v＜k，v＝v+1，跳转至5.2步；否则，执行5.11步；

5.11如果u＜M，u＝u+1，跳转至5.2步；否则，执行第六步。

第六步，如果j＝j+1，跳转至第五步。否则，j＝j+1，执行第七步。

第七步，如果

i＝i+1，跳转至第五步。否则，结束。

如果SIMD部件包含乘加指令，5.6和5.7可以合并为一个步骤。

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

采用本发明可以实现任意维度(第一个输入矩阵A的列数等于第二个输入矩阵B的行数)的两个矩阵在双精度SIMD部件上相乘。本发明中5.3步对A_i是按行进行访问的(5.10步增加v，其后的5.11步增加u，因此5.3步的a_uv是按行访问的)，5.5步中对B_j也是按行访问的，即对输入矩阵A和B可以按照其在存储中的相同顺序访问(如果A和B是按列存储的，将5.10和5.11调换顺序，5.5中取B_j每列的W个数据，就可以实现对输入矩阵A和B都是进行按列访问)，避免了对其中的一个矩阵进行转置操作。同时，在子矩阵计算过程中，矩阵B的内容可以被重复利用，减少了访问矩阵B的时间，提高了双精度SIMD部件的计算效率。

附图说明

图1是本发明总体流程图；

图2是传统的面向SIMD部件的矩阵乘方法的一个实例。

图3是使用本发明的方法进行矩阵乘的一个实例。

具体实现方式

图1是本发明总体流程图，本发明的总体过程为：

第一步，增加矩阵B的列数并进行数据填充；

第二步，增加矩阵C的列数，并将C的内容初始化为全0；

第三步，矩阵B进行分块；

第四步，矩阵A进行分块；

第五步，子矩阵相乘，并将结果加到矩阵C中的行的结果和上；

第六步，是否遍历了B的位于同行的分块，如果是，执行第七步；否则，跳转至第五步；

第七步，是否遍历了A的所有分块，如果是，程序结束；否则，跳转至第五步。

为了检验面向SIMD部件的矩阵乘实现效果，使用国防科大飞腾CPU为实现平台，飞腾CPU的SIMD部件的向量长度为4，向量寄存器的个数为32。在该平台上采用C语言实现了SIMD部件的矩阵乘方法。假设两个输入矩阵是64×64的矩阵，根据本发明，B被划分为64个子块B_j，每个子块的大小为16×4；A被划分为4个子块A_i，每个子块的大小为64×16。图3给出了使用本发明实现两个64×64矩阵乘的方法，子块A_i(i为从1至4的整数)分别和子块B_j(j为从1+16*(i-1)至16*i的整数)进行相乘，将子矩阵相乘结果进行累加得到结果矩阵。因此采用本发明进行A和B相乘是对A和B按相同顺序进行访问，无需对A或B进行转置操作。图2给出了面向SIMD部件的两个64×64矩阵乘的传统方法，需要对B进行转置后进行矩阵乘的计算。

在飞腾CPU上，使用传统方法，转置B矩阵的时间为0.002秒，SIMD部件的计算时间为0.056秒，矩阵乘总共的计算时间为0.058秒。使用本发明时，面向SIMD部件的矩阵乘计算时间为0.055秒，性能提升了5.2％。

同时，对640×640的两个矩阵乘进行了实验，使用传统方法，转置B矩阵的时间为0.033秒，SIMD部件的计算时间为0.82秒，矩阵乘总共的计算时间为0.853秒；使用本发明时，矩阵乘计算时间为0.81秒，性能提升了5.04％。

Claims (1)

1.一种面向双精度SIMD部件的矩阵乘实现方法，其特征在于包括以下步骤：

第一步，对于输入矩阵A和B，当N不是W的整数倍时，增加输入矩阵B的列数，将B增加W-N％W列，使B的列数为N+W-N％W，％表示模运算，增加的列使用0进行数据填充；A矩阵的大小是M*K，B矩阵为K*N，A、B均为双精度浮点矩阵，W为双精度SIMD部件的向量长度，即一个向量包含W个双精度浮点数；N、W、K全是整数；结果矩阵为C，大小为M*N；

第二步，增加结果矩阵C的列数，使C的列数为N+W-N％W，并将矩阵C的内容初始化为全0；

第三步，根据SIMD部件的向量寄存器数目VN对矩阵B进行分块，将K*N的矩阵B划分成k*n的子块Bj，

其中n必须是W的整数倍；VN为正整数；

子矩阵B_j的大小即n和k值须满足：

2+n*(k+1)/W＜VN且n％W＝0和k％W＝0，

2+n*(k+1)/W为每次子矩阵相乘所需的最少的向量寄存器数目；

第四步，将M*K的矩阵A划分成M*k的子块A_i，

表示下取整；

第五步，子矩阵A_i和B_j在SIMD部件中相乘，并将结果累加到结果矩阵C中；

5.1令i＝1，j＝1，u＝1，v＝1；

5.2将n/W个结果向量寄存器V_s的内容初始化为0，1≤s≤n/W；

5.3使用广播指令ld1toW将A_i中的一个元素a_uv取至向量寄存器V₀；

5.4令P＝1；

5.5如果u等于1，使用向量访存指令将B_j中的第v行元素中从(P-1)*W+1至P*W的元素取至向量寄存器V_Z中，1+n/W≤z≤n*(k+1)/W，执行第5.6步；如果u不等于1，则数据已经存放在向量寄存器中，执行第5.6步；

5.6V₀和V_Z进行向量乘法操作，将结果存储在向量寄存器V_t中t＝1+n*(k+1)/W；

5.7V_t和结果向量寄存器V_s进行向量加操作，将结果存放在V_s中；

5.8如果P＜n/W，P＝P+1，跳转至5.5；否则，执行5.9步；

5.9将n/W个结果向量寄存器V_s中的数据和C的第u行中第(u-1)*n+1+(i-1)*n列至第u*n+(i-1)*n列的n个数据进行累加，并将结果写到C中；

5.10如果v＜k，v＝v+1，跳转至5.2步；否则，执行5.11步；

5.11如果u＜M，u＝u+1，跳转至5.2步；否则，执行第六步；

第六步，如果j＝j+1，跳转至第五步，否则，j＝j+1，执行第七步；

第七步，如果

i＝i+1，跳转至第五步，否则，结束。

Images