并行计算.5矩阵.ppt

《并行计算.5矩阵.ppt》由会员分享，可在线阅读，更多相关《并行计算.5矩阵.ppt（78页珍藏版）》请在三一文库上搜索。

1、并 行 计 算,稠密矩阵运算,稠密矩阵运算,矩阵是一种重要的数据结构，它可以表达很丰富的内容。矩阵可以分为两类，一类是稠密（dense）矩阵，其特点是多数元素都非零；第二类矩阵称为稀疏矩阵（sparse matrices），特点是多数元素都为零。之所以这样区分，是因为稀疏矩阵可以通过只存储非零元素（及其索引）而节省空间，并且稀疏矩阵在很多操作上可以优化。,矩阵的划分,为了对一个矩阵进行并行操作，必须将其进行划分并映射到不同的处理器上，尤其是科学计算中矩阵的规模一般都很庞大 的情况下更是如此。以后将会看到，不同的矩阵划分方式对算法的性能有着不容忽视的影响，因此确定那些划分和映射方式最有利于性能提

2、高是很重要的。在这儿将讨论常用的矩阵划分和映射方式。,带状划分,带状划分就是把矩阵按照行或列分成几部分，分别映射到各个处理器。如果分到每个处理器的各行或列是连续的，则称为块带状划分（Block-Striped）；相对的，如果是按照行号或者列号取模而进行的矩阵划分则称为循环带状划分（Cyclic-Striped）。 下图是一个1616的矩阵带状划分到4各处理器的例子，左右分别为列方向的块带状划分和行方向的循环带状划分。 带状划分最多能够把一个nn的矩阵划分到n的处理器上。,带状划分,1616阶矩阵，p=4 列带状划分 行循环带状划分,带状划分,3，27 27矩阵的3种带状划分,棋盘划分,棋盘划分

3、把矩阵划分称为若干个子矩阵并映射到不同的处理器上，每个处理器都不包含完整的行和列。和条带状划分类似，棋盘划分也可以分为块棋盘划分和循环棋盘划分 对于实际的网格状互联网络结构，矩阵的棋盘划分可以方便地映射到其上。比如把一个nn的矩阵按照棋盘划分映射到一个qq的二维网格上（设p=qq），则每个处理器包含 个矩阵元素。 利用棋盘划分可以划分到更细的粒度（最多可以划分成nn块），从而可以达到更高的并行度。,棋盘划分,88阶矩阵，p=16 块棋盘划分 循环棋盘划分,棋盘划分,示例：p4，1616矩阵的3种棋盘划分,矩阵的转置,转置（Transposition）是基本的矩阵运算。一个矩阵A的转置记为 ，它

4、是将矩阵A的元素沿对角线互换而得到的。矩阵转置的串行算法很简单，只需要把上三角（不包括对角线）的元素循环一遍，每个元素与其对称位置的元素交换位置即可，整个过程只需要一个单位的多余空间，时间复杂度为O(n2)。下面讨论在不同的矩阵划分方式下的并行矩阵转置算法。,带状划分的矩阵转置,划分: Ann分成p个(n/p)n大小的带 算法: Pi有p-1个(n/p)(n/p)大小子块发送到另外p-1个处理器中; 每个处理器本地交换相应的元素,棋盘划分的矩阵转置,仅讨论网格块棋盘划分的矩阵转置算法，循环棋盘划分和映射同样考虑即可。下面分两种不同的网络互联结构分别讨论：二维网格互联结构、超立方体结构。,网孔连

5、接,划分: 算法: 按mesh连接进行块转置(不同处理器间) 进行块内转置(同一处理器内),网孔连接,情形1: p=n2。 通讯步 转置后,44矩阵转置,44矩阵转置,网孔连接,设处理器数目pn2，则整个nn的矩阵分成p个大小为 的子块，这些子块自然影射到网格互联结构的p个处理器上。算法可以分两步来完成：第一步，将各个子块看作一个整体进行图 (a)所示的子块转置，这一步由各处理器互相通信；第二步，各处理器分别对自己的内部子块进行局部转置，如图 (b)所示。,情形2:,网孔连接,通讯步 转置后,网孔连接,整个算法的运行时间分析如下：第一步子块间转置最长的路径约为 ，移动单个子块的时间为 ，所以各

6、个子块到达目的地的时间为 ；第二步，大小为 的子块转置需要的时间量级是n2/p，所以算法的总运行时间为：,超立方连接,转置也可以递归地进行，称为递归转置算法（Recursive Transposition Algorithm）。基本思路是一个nn的矩阵先分成4个子矩阵，把每个子矩阵看作一个整体对整个大矩阵进行转置；然后每个子矩阵再分成4个更小的矩阵，逐次递归,超立方连接,这种递归跟超立方体的递归构建过程很像。利用这一点可以在超立方体上实现矩阵转置。先把nn的矩阵分成4个(n/2)(n/2) 的子矩阵，相应地，一个具有p个处理器的超立方体可以看作是4个p/4超立方体所构成。递归一直到每个子超立方

7、体只含有一个处理器为止。下面来分析算法的运行 时间，经过 次后递归结束，此时各个子块的大小为 ，这些子块内部转置的时间量级为n2/p。每个待通信的子块大小为n2/p，使用存储转发选路所需的时间为2( )，递归 次的总通信时间为 。这样，算法总的运行时间为：,超立方连接,划分: 算法: 对Aij递归应用进行转置，直至分块矩阵的元素处于同一处理器； 进行同一处理器的内部转置。 运行时间:,超立方连接：示例,矩阵向量乘法,一个矩阵与向量相乘，结果得到一个向量。矩阵向量相乘的串行算法具有O(n2)的时间复杂度。不同的矩阵划分与映射方式，以及不同的网络互联结构所产生的并行化效果也不同。,带状划分的矩阵-

8、向量乘法,=n时，最初每个处理器存储矩阵的一行以及向量的一个元素。由于矩阵与向量相乘需要用到整个向量，所以每个处理器都要把自己所存储的向量元素发送给其它处理器，这是一个多到多的通信过程。向量分布到各个处理器上后就可以进行乘法运算了。,带状划分的矩阵-向量乘法,处理器数目pn的情况类似，只是每个处理器存放矩阵的n/p行和向量的n/p个元素。同样也很容易理解，块带状划分与循环带状划分是没有本质区别的。,带状划分的矩阵-向量乘法,划分(行带状划分): Pi存放xi和ai,0,ai,1,ai,n-1, 并输出yi 算法: 对p=n情形 每个Pi向其他处理器播送xi(多到多播送)； 每个Pi计算； 注:

9、 对pn情形,算法中Pi要播送X中相应的n/p个分量 (1)超立方连接的计算时间 (2)网孔连接的计算时间,带状划分的矩阵-向量乘法,示例,带状划分的矩阵-向量乘法,初始划分,带状划分的矩阵-向量乘法,多到多通信传输向量信息,带状划分的矩阵-向量乘法,多到多通信之后的向量分布,带状划分的矩阵-向量乘法,最终结果,棋盘划分的矩阵-向量乘法,只讨论块棋盘划分和映射的情况，循环棋盘划分和映射的情况类似。 从最极端的情况p=n2开始，此时每个处理器存放矩阵的一个元素，最初某些（n个）处理器中存放有向量的一个元素。设向量X的元素Xi存放在处理器Pii中，如果不是这样，则本行的某处理器通过一对一通信将向量

10、元素传送到Pii。然后处理器Pii需要把元素Xi传送给本列的其它处理器，这样每个处理器就都可以进行并行地进行元素相乘了。元素相乘的结果再按照行累加起来，并交给每行的一个处理器。,棋盘划分的矩阵-向量乘法,划分(块棋盘划分): Pij存放ai,j, xi置入Pi,i中 算法: 对p=n2情形 每个Pi,i向Pj,i播送xi(一到多播送)； 按行方向进行乘-加与积累运算，最后一列Pi,n-1收集的结果为yi； 注: 对pn2情形,p个处理器排成 的二维网孔， 算法中Pi,i向Pj,i播送X中相应的 个分量 网孔连接的计算时间Tp(CT): .X中相应分量置入Pi,i的通讯时间: .按列一到多播送时

11、间: .按行单点积累的时间:,棋盘划分的矩阵-向量乘法,示例,棋盘划分的矩阵-向量乘法,初始划分和通信,棋盘划分的矩阵-向量乘法,沿着列方向传送元素,棋盘划分的矩阵-向量乘法,元素相乘和结果按行累加,棋盘划分的矩阵-向量乘法,结果向量的分布,带状与棋盘划分比较,以网孔链接为例 网孔上带状划分的运行时间 网孔上棋盘划分的运行时间 棋盘划分要比带状划分快。,矩阵乘法,在开始并行化之前，先引入分块矩阵操作的概念。在进行矩阵操作时，我们总是倾向于把矩阵操作转化为对矩阵元素的操作。其实，可以把一个子矩阵块看作一个超 元素，把矩阵操作转化为对矩阵块（超元素）的操作；如果有必要，矩阵块还可以进一步分解为更小

12、的矩阵块。这种概念叫做分块矩阵操作。比如一个nn的矩阵 可以看作由qq个矩阵块所构成，每个矩阵块是一个(n/q)(n/q)的矩阵。,矩阵乘法符号及定义,A中元素的第1下标与B中元素的第2下标相一致（对准）,矩阵乘法并行实现方法,计算结构：二维阵列 空间对准(元素已加载到阵列中） Cannons , Foxs，DNS 时间对准(元素未加载到阵列中） Systolic,简单并行分块乘法,分块: A、B和C分成 的方块阵Ai,j、Bi,j和Ci,j, 大小均为 p个处理器编号为 , Pi,j存放Ai,j、Bi,j和Ci,j。 算法: 通讯:每行处理器进行A矩阵块的多到多播送(得到Ai,k, k=0

13、) 每列处理器进行B矩阵块的多到多播送(得到Bk,j, k=0 ) 乘-加运算: Pi,j做 运行时间 (1)超立方连接： 的时间 的时间,简单并行分块乘法,运行时间 (2)二维环绕网孔连接： 的时间： 的时间t2=n3/p 注 (1)本算法的缺点是对处理器的存储要求过大 每个处理器有 个块，每块大小为n2/p, 所以需要 , p个处理器共需要 ， 是串行算法的 倍 (2)p=n2时，t(n)=O(n), c(n)=O(n3),Cannon乘法,Cannon乘法的矩阵分块方式与简单分块矩阵法完全一样，但它没有等到所有的块Aik与Bkj（ ）都存储到本机后才进行Cij的计算，而是每拿到一个Aik

14、和Bkj就进行乘法计算，小块计算完毕后就可以腾出空间给别的矩阵块使用。,Cannon乘法,算法的分为准备阶段和循环阶段，在准备阶段，块Aij向左循环移动i步，Bij向上循环移动j步。在循环阶段的每个循环中，处理器Pij计算它上面两个当前块Aij和Bij的乘积，然后将Aij向左循环移动一步，Bij向上循环移动一步。所谓循环移动一步，就是移动一步之后取模，与算术移动相对应（可以联想到汇编语言中的循环移位）。,Cannon乘法,分块: A、B和C分成 的方块阵Ai,j、Bi,j和Ci,j, 大小 均为p个处理器编号为 , Pi,j存放Ai,j、Bi,j和Ci,j(n p),Cannon乘法,算法原理

15、 (非形式描述) 所有块Ai,j(0i,j )向左循环移动i步(按行移位)； 所有块Bi,j(0i,j )向上循环移动j步(按列移位)； 所有处理器Pi,j做执行Ai,j和Bi,j的乘-加运算; A的每个块向左循环移动一步; B的每个块向上循环移动一步; 转执行 次;,Cannon乘法,示例： A44， B44, p=16,Cannon乘法,示例： A44， B44, p=16,Cannon乘法,示例： A44， B44, p=16,Cannon乘法,上图演示了在16个处理器上完成两个44矩阵A和B相乘的过程。在准备阶段，矩阵A的第0行（A00到A03）不动，即向左移动0步； A10 至A13

16、 向左循环移动一位； A20 至A23 向左循环移动2位； A30 至A33 向左循环移动3位。矩阵B的第0列（B00到B30）向上移动0位；B01至B31向上循环移动一位；B02至B32向上循环移动2位；B03至B33向上循环移动3位。在循环阶段，每计算一次块矩阵的乘法就把用过的块矩阵A移动给左邻居，把块矩阵B移动给上邻居。,Cannon乘法,算法描述: Cannon分块乘法算法 /输入: Ann, Bnn; 输出: Cnn Begin (1)for k=0 to do for all Pi,j par-do (i) if ik then Ai,j Ai,(j+1)mod endif (ii

17、)if jk then Bi,j B(i+1)mod , j endif endfor endfor (2)for all Pi,j par-do Ci,j=0 endfor,(3)for k=0 to do for all Pi,j par-do (i) Ci,j=Ci,j+Ai,jBi,j (ii) Ai,j Ai,(j+1)mod (iii)Bi,j B(i+1)mod , j endfor endfor End,时间分析：,Fox乘法,Fox算法同样通过循环移位的办法来达到节省存储空间的目的，但与Cannon算法不同的是，它通过一对多广播的方式处理A的子矩阵块，而象 Cannon那样的

18、用循环移位方法处理B的子矩阵块。根据对称性，也可以对矩阵A的子阵进行循环移位而对矩阵B的子阵采取广播的方式，效果一样。设处理器个 数 ，则算法的要点如下：,Fox乘法,1. 所选中的对角线Aii向所在行的q个处理器进行一对多广播； 2. 各处理器将自己所拥有的A和B的子块进行矩阵相乘运算； 3. B矩阵的块向上循环移动一位，从下面接受一个新的B矩阵块； 4. 选择A的一个矩阵块作为广播源，选择方法是：如果Aij是上次的广播源，则本次的广播源是 。其中%表示取模运算。转步骤1。,Fox乘法,分块: 同Cannon分块算法 算法原理 Ai,i向所在行的其他处理器 进行一到多播送； 各处理器将收到的

19、A块与原 有的B块进行乘-加运算; B块向上循环移动一步; 如果Ai,j是上次第i行播送的块， 本次选择 向所在行的其他处理器进行一到多播送; 转执行 次;,Fox乘法,示例： A44， B44, p=16 (a) (b),Fox乘法,示例： A44， B44, p=16 (c) (d),B2,0,B3,0,B2,1,B3,1,B0,1,B3,2,B0,2,B1,2,B2,3,B0,3,B1,3,B3,0,B1,0,B3,1,B0,1,B2,1,B3,2,B1,2,B0,3,B2,3,B0,2,B1,3,B2,0,A0,2 B2,2,A1,3 B3,3,A2,0 B0,0,B1,0,A3,1

20、B1,1,A0,3 B3,3,A1,0 B0,2,A2,1 B1,1,A3,2 B2,2,Systolic乘法,Systolic乘法,a1,4,b4,1,b3,1,b2,1,b2,2,b4,2,b3,2,b2,3,b3,3,b4,3,b2,4,b3,4,b4,4,a1,3,a1,1,a1,2,a2,4,a2,1,a2,2,a2,3,a3,1,a3,2,a3,3,a3,4,b1,1,b1,2,b1,3,b1,4,Step 1,P1,1 c1,1,P1,2 c1,2,P1,3 c1,3,P1,4 c1,4,P2, 1 c2,1,P2,2 c2,2,P2,3 c2,3,P2,4 c2,4,P3,1

21、c3,1,P3,2 c3,2,P3,3 c3,3,P3,4 c3,4,Systolic乘法,c1,1,c1,2,c1,3,c1,4,c2,1,c2,2,c2,3,c2,4,c3,1,c3,2,c3,3,c3,4,b3,1,b2,1,b2,2,b4,2,b3,2,b2,3,b3,3,b4,3,b2,4,b3,4,b4,4,a1,3,a1,1,a1,2,a2,4,a2,1,a2,2,a2,3,a3,1,a3,2,a3,3,a3,4,b1,1,b1,2,b1,3,b1,4,a1,4,b4,1,+,Step 2,Systolic乘法,c1,1,c1,2,c1,3,c1,4,c2,1,c2,2,c2,3

22、,c2,4,c3,1,c3,2,c3,3,c3,4,b2,1,b2,2,b3,2,b2,3,b3,3,b4,3,b2,4,b3,4,b4,4,a1,1,a1,2,a2,1,a2,2,a2,3,a3,1,a3,2,a3,3,a3,4,b1,1,b1,2,b1,3,b1,4,a1,3,b3,1,+,a1,4,b4,2,+,a2,4,b4,1,+,Step 3,Systolic乘法,c1,1,c1,2,c1,3,c1,4,c2,1,c2,2,c2,3,c2,4,c3,1,c3,2,c3,3,c3,4,b2,2,b2,3,b3,3,b2,4,b3,4,b4,4,a1,1,a2,1,a2,2,a3,1,

23、a3,2,a3,3,b1,1,b1,2,b1,3,b1,4,a1,2,b2,1,+,a1,3,b3,2,+,a2,3,b3,1,+,a1,4,b4,3,+,a3,4,b4,1,+,a2,4,b4,2,+,Step 4,Systolic乘法,c1,1,c1,2,c1,3,c1,4,c2,1,c2,2,c2,3,c2,4,c3,1,c3,2,c3,3,c3,4,b2,3,b2,4,b3,4,a2,1,a3,1,a3,2,b1,2,b1,3,b1,4,a1,1,b1,1,+,a1,2,b2,2,+,a2,2,b2,1,+,a1,3,b3,3,+,a3,3,b3,1,+,a2,3,b3,2,+,a1,

24、4,b4,4,+,a2,4,b4,3,+,a3,4,b4,2,+,Step 5,Systolic乘法,c1,1,c1,2,c1,3,c1,4,c2,1,c2,2,c2,3,c2,4,c3,1,c3,2,c3,3,c3,4,b2,4,a3,1,b1,3,b1,4,a1,1,b1,2,+,a2,1,b1,1,+,a1,2,b2,3,+,a3,2,b2,1,+,a2,2,b2,2,+,a1,3,b3,4,+,a2,3,b3,3,+,a3,3,b3,2,+,a2,4,b4,4,+,a3,4,b4,3,+,Step 6,Systolic乘法,c1,1,c1,2,c1,3,c1,4,c2,1,c2,2,c

25、2,3,c2,4,c3,1,c3,2,c3,3,c3,4,b1,4,a1,1,b1,3,+,a3,1,b1,1,+,a2,1,b1,2,+,a1,2,b2,4,+,a2,2,b2,3,+,a3,2,b3,2,+,a2,3,b3,4,+,a3,3,b3,3,+,a3,4,b4,4,+,Step 7,Systolic乘法,c1,1,c1,2,c1,3,c1,4,c2,1,c2,2,c2,3,c2,4,c3,1,c3,2,c3,3,c3,4,a1,1,b1,4,+,a2,1,b1,3,+,a3,1,b1,2,+,a2,2,b2,4,+,a3,2,b2,3,+,a3,3,b3,4,+,Step 8,S

26、ystolic乘法,c1,1,c1,2,c1,3,c1,4,c2,1,c2,2,c2,3,c2,4,c3,1,c3,2,c3,3,c3,4,a2,1,b1,4,+,a3,1,b1,3,+,a3,2,b2,4,+,Step 9,Systolic乘法,c1,1,c1,2,c1,3,c1,4,c2,1,c2,2,c2,3,c2,4,c3,1,c3,2,c3,3,c3,4,a3,1,b1,4,+,Step 10,Systolic乘法,P1,1 c1,1,P1,2 c1,2,P1,3 c1,3,P1,4 c1,4,P2, 1 c2,1,P2,2 c2,2,P2,3 c2,3,P2,4 c2,4,P3,1

27、 c3,1,P3,2 c3,2,P3,3 c3,3,P3,4 c3,4,Over,c1,1 = a1,1 b1,1 + a1,2 b2,1 + a1,3 b3,1 + a1,4 b4,1 c1,2 = a1,1 b1,2 + a1,2 b2,2 + a1,3 b3,2 + a1,4 b4,2 c3,4 = a3,1 b1,4 + a3,2 b2,4 + a3,3 b3,4 + a3,4 b4,4,Systolic乘法,Systolic算法 /输入: Amn, Bnk; 输出: Cmk Begin for i=1 to m par- do for j=1 to k par-do (i) ci,j

28、 = 0 (ii) while Pi,j 收到a和b时 do ci,j = ci,j +ab if i m then 发送b给Pi+1,j endif if j k then 发送a给Pi,j+1 endif endwhile endfor endfor End,DNS乘法,背景: 由Dekel、Nassimi和Sahni提出的SIMD-CC上的矩阵乘法, 处理器 数目为n3, 运行时间为O(logn), 是一种速度很快的算法。 基本思想: 通过一到一和一到多的播送办法，使得处理器(k,i,j)拥有ai,k和bk,j, 进行本地相乘,再沿k方向进行单点积累求和,结果存储在处理器(0,i,j)中

29、。 处理器编号: 处理器数p=n3= (2q)3=23q, 处理器Pr位于位置(k,i,j), 这里r=kn2+in+j, (0i, j, kn-1)。位于(k,i,j)的处理器Pr的三个寄存器 Ar,Br,Cr分别表示为Ak,i,j, Bk,i,j和Ck,i,j, 初始时均为0。 算法: 初始时ai,j和bi,j存储于寄存器A0,i,j和B0,i,j; 数据复制:A,B同时在k维复制(一到一播送)； A在j维复制(一到多播送); B在i维复制(一到多播送); 相乘运算:所有处理器的A、B寄存器两两相乘; 求和运算:沿k方向进行单点积累求和;,示例 C00=1(-5)+27=9 C01=1(-

30、6)+28=10 C10=3(-5)+47=13 C11=3(-6)+48=14,初始加载,(b)A,B沿k维复制,(c)A沿j维复制,(d)B沿i维复制,(e)点积,(f)沿k维求和,B,B,A,A,000,010,001,011,100,110,101,111,P0,P1,P2,P3,P4,P5,P6,P7,算法描述: /令r(m)表示r的第m位取反； /p, rm=d表示r(0rp-1)的集合, 这里r的二 /进制第m位为d; /输入: Ann, Bnn; 输出: Cnn Begin /以n=2, p=8=23举例, q=1, r=(r2r1r0)2 (1)for m=3q-1 to 2

31、q do /按k维复制A,B, m=2 for all r in p, rm=0 par-do /r2=0的r (1.1) Ar(m) Ar /A(100)A(000)等 (1.2) Br(m) Br /B(100)B(000)等 endfor endfor (2)for m=q-1 to 0 do /按j维复制A, m=0 for all r in p, rm= r2q+m par-do /r0=r2的r Ar(m) Ar /A(001)A(000),A(100)A(101) endfor /A(011)A(010),A(110)A(111) endfor,(3)for m=2q-1 to

32、q do /按i维复制B,m=1 for all r in p, rm= rq+m par-do/r1=r2的r Br(m) Br /B(010)B(000),B(100)B(110) endfor /B(011)B(001),B(101)B(111) endfor (4)for r=0 to p-1 par-do /相乘, all Pr Cr=ArBr endfor (5)for m=2q to 3q-1 do /求和,m=2 for r=0 to p-1 par-do Cr=Cr+Cr(m) endfor endfor End,DNS乘法,示例,DNS乘法,DNS乘法,DNS乘法,图中，把

33、处理器分成n个平面，每个平面上有nn个处理器，从图上可以看到，每个平面对应一个不同的k值。 处理器Pr计算行AI, *和列B*, j的点积，计算所需要的数据通过处理器之间的通信传送到目标处理器，使得处理器Pk,i,j上存储有数据Ai, k和Bk, j。 把数据放置到合适的处理器上可以通过两个通信步骤来实现。第一步把数据从k=0的平面复制到其它平面；第二步各个平面内部的数据扩散。第一步的数据复制机制可以用一句话来概括：把A*, j复制到第k=j的平面；把Bi, *复制到k=i的平面。第二步的数据扩散方式是：矩阵A的元素沿着j方向一对多通信复制；矩阵B的元素沿着i方向一对多通信复制。 经过这两个通信步骤后，Ai, k和Bk, j的乘积就可以在处理器Pk,i,j上进行了。处理器P*,i,j上的乘积求和后在P0,i,j上收集。,