《并行计算》课程教学资源（讲义）矩阵运算

14 Cannon乘法 141 Cannon乘法的原理 Cannon算法是一种存储有效的算法。为了使两矩阵下标满足相乘的要求,它和上一节 的并行分块乘法不同,不是仅仅让B矩阵的各列块循环移动,而是有目的地让A的各行块 以及B的各列块皆施行循环移位,从而实现对C的子块的计算。将矩阵A和B分成p个方 块A和Bn,(0≤,j≤VP-1,每块大小为「mx「mF,并将它们分配给√Px√P个处 理器(P00,Po1,P p-1√p-1 开始时处理器Pg存放块A和B,并负责计算块Cy,然 后算法开始执行: (1)将块A向左循环移动i步;将块B向上循环移动j步; (2)P执行乘加运算后将块A向左循环移动1步,块Bn向上循环移动1步 (3)重复第(2)步,总共执行√P次乘加运算和√P次块A1和Bn的循环单步移位 142 Cannon乘法的并行算法 图13示例了在16个处理器上,用 Cannon算法执行A×x4XB4x4的过程。其中(a)和(b) 对应于上述算法的第(1)步;(c)、(d)、(e)、(f)对应于上述算法的第(2)和第(3)步。在算法第(l) 步时,A矩阵的第0列不移位,第Ⅰ行循环左移1位,第2行循环左移2位,第3行循环左 移3位:类似地,B矩阵的第0行不移位,第1列循环上移Ⅰ位,第2列循环上移2列,第 3列循环上移3列。这样 Cannon算法具体描述如下: 算法187 Cannon乘法算法 输入:Anxn,Bn 输出:C 对所有处理器 my_ rank( my rank=0,…p-1)同时执行如下的算法 (1)计算子块的行号r= my rank/sqrt(p) 计算子块的列号产 my rank mod sqrt(p) if(p>k) then Leftmoveonestep(a) end if/*a循环左移至同行相邻处理器中* if(>k)then Up tep(b) end if/b循环上移至同列相邻处理器中* end for (4 )for k=0 to p for i=0 to m-1 d

1.4 Cannon 乘法 1.4.1 Cannon 乘法的原理 Cannon 算法是一种存储有效的算法。为了使两矩阵下标满足相乘的要求，它和上一节 的并行分块乘法不同，不是仅仅让 B 矩阵的各列块循环移动，而是有目的地让 A 的各行块 以及 B 的各列块皆施行循环移位，从而实现对 C 的子块的计算。将矩阵 A 和 B 分成 p 个方 块 Aij和 Bij，(0  i, j  p −1) ，每块大小为 n/ p   n/ p  ，并将它们分配给 p  p 个处 理器 ( , ,..., ) 00 01 p−1 p−1 P P P 。开始时处理器 Pij 存放块 Aij 和 Bij，并负责计算块 Cij，然 后算法开始执行： ⑴将块 Aij 向左循环移动 i 步；将块 Bij 向上循环移动 j 步； ⑵Pij 执行乘加运算后将块 Aij 向左循环移动 1 步，块 Bij 向上循环移动 1 步； ⑶重复第⑵步，总共执行 p 次乘加运算和 p 次块 Aij 和 Bij 的循环单步移位。 1.4.2 Cannon 乘法的并行算法 图 1.3 示例了在 16 个处理器上，用 Cannon 算法执行 A4×4×B4×4 的过程。其中(a)和(b) 对应于上述算法的第⑴步；(c)、(d)、(e)、(f)对应于上述算法的第⑵和第⑶步。在算法第⑴ 步时，A 矩阵的第 0 列不移位，第 1 行循环左移 1 位，第 2 行循环左移 2 位，第 3 行循环左 移 3 位；类似地，B 矩阵的第 0 行不移位，第 1 列循环上移 1 位，第 2 列循环上移 2 列，第 3 列循环上移 3 列。这样 Cannon 算法具体描述如下： 算法 18.7 Cannon 乘法算法 输入：An×n，Bn×n 输出：Cn×n Begin 对所有处理器 my_rank(my_rank=0,…,p-1)同时执行如下的算法: (1)计算子块的行号 i=my_rank/sqrt(p) 计算子块的列号 j=my_rank mod sqrt(p) (2)for k=0 to p -1 do if (i>k) then Leftmoveonestep(a) end if /* a 循环左移至同行相邻处理器中*/ if (j>k) then Upmoveonestep(b) end if /* b 循环上移至同列相邻处理器中*/ end for (3)for i=0 to m-1 do for j=0 to m-1 do c[i,j]=0 end for end for (4)for k=0 to p -1 do for i=0 to m-1 do for j=0 to m-1 do for k1=0 to m-1 do

c[∥=c[)+a[k1]*bA1 Leftmoveonestep(a)/子块a循环左移至同行相邻的处理器中* Upmoveonestep(b)/*子块b循环上移至同列相邻的处理器中* end for End A3 (a)阵起始对准 b)B阵起始对准 (c)对准后的A和B (d)第一次移位后的子阵位置 B3 303 ::2:2:签2:2 (e)第二次移位后的子阵位置 (f)第三次移位后的子阵位置 图1.316个处理器上的 Cannon乘法过程 这里函数 Leftmoveonestep(a)表示子块a在编号处于同行的处理器之间以循环左移的方 式移至相邻的处理器中;函数 Upmoveonestep(b)表示子块b在编号处于同列的处理器之间以 循环上移的方式移至相邻的处理器中。这里我们以函数 Leftmoveonestep(a)为例,给出处理 器间交换数据的过程:

c[i,j]= c[i,j]+ a[i,k1]* b[k1,j] end for end for end for Leftmoveonestep(a) /*子块 a 循环左移至同行相邻的处理器中*/ Upmoveonestep(b) /*子块 b 循环上移至同列相邻的处理器中*/ end for End (a)A阵起始对准 (b)B阵起始对准 0,0 1,0 2,0 3,0 0,1 1,1 2,1 3,1 0,2 1,2 2,2 3,2 0,3 1,3 2,3 3,3 0,0 1,0 2,0 3,0 0,1 1,1 2,1 3,1 0,2 1,2 2,2 3,2 0,3 1,3 2,3 3,3 A A A A B B B B A A A A B B B B A A A A B B B B A A A A B B B B 0,1 A 0,0 A 0,3 A 0,2 A 1,2 A 1,1 A 1,0 A 1,3 A 2,3 A 2,2 A 2,1 A 2,0 A 3,0 A 3,3 A 3,2 A 3,1 A 0,0 B 1,0 B 2,0 B 3,0 B 0,1 B 1,1 B 2,1 B 3,1 B 0,2 B 1,2 B 2,2 B 3,2 B 0,3 B 1,3 B 2,3 B 3,3 B 0,2 A 0,1 A 0,0 A 0,3 A 1,3 A 1,2 A 1,1 A 1,0 A 2,0 A 2,3 A 2,2 A 2,1 A 3,1 A 3,0 A 3,3 A 3,2 A 1,0 B 2,0 B 3,0 B 0,0 B 1,1 B 2,1 B 3,1 B 0,1 B 1,2 B 2,2 B 3,2 B 0,2 B 1,3 B 2,3 B 3,3 B 0,3 B (c)对准后的A和B (d)第一次移位后的子阵位置 0,3 A 0,2 A 0,1 A 0,0 A 1,0 A 1,3 A 1,2 A 1,1 A 2,1 A 2,0 A 2,3 A 2,2 A 3,2 A 3,1 A 3,0 A 3,3 A 2,0 B 3,0 B 0,0 B 1,0 B 2,1 B 3,1 B 0,1 B 1,1 B 2,2 B 3,2 B 0,2 B 1,2 B 2,3 B 3,3 B 0,3 B 1,3 B (e)第二次移位后的子阵位置 0,3 A 1,0 A 2,1 A 3,2 A (f)第三次移位后的子阵位置 0,0 B 1,0 B 2,0 B 3,0 B 0,1 A 0,0 A 0,2 A 1,2 A 1,1 A 1,3 A 2,3 A 2,2 A 2,0 A 3,0 A 3,3 A 3,1 A 0,1 B 1,1 B 2,1 B 3,1 B 0,2 B 1,2 B 2,2 B 3,2 B 0,3 B 1,3 B 2,3 B 3,3 B 图 1.3 16 个处理器上的 Cannon 乘法过程 这里函数 Leftmoveonestep(a)表示子块 a 在编号处于同行的处理器之间以循环左移的方 式移至相邻的处理器中；函数 Upmoveonestep(b)表示子块 b 在编号处于同列的处理器之间以 循环上移的方式移至相邻的处理器中。这里我们以函数 Leftmoveonestep(a)为例，给出处理 器间交换数据的过程：

算法18.8函数 Leftmoveonestep(a)的基本算法 (1)if(=0)then/*最左端的子块* (1.1)将所存的A的子块发送到同行最右端子块所在的处理器中 (1.2)接收其右邻处理器中发来的A的子块 end if (2)if(G=sqrt(p}-1)and(mod2=0)then/最右端子块处理器且块列号为偶数* (2.1)将所存的A的子块发送到其左邻处理器中 (2.2)接收其同行最左端子块所在的处理器发来的A的子块 end if (3)if(G=sqrt(p}-1)and(mod2≠0)then/最右端子块处理器且块列号为奇数* (3,1)将所存的A的子块在缓冲区buer中做备份 (3,2)接收其同行最左端子块所在的处理器发来的A的子块 (33)将在缓冲区 buffer中所存的A的子块发送到其左邻处理器中 end if (4)if(≠sqrt(p)-)and(md2=0)and(≠0)then/其余的偶数号处理器 (41)将所存的A的子块发送到其左邻处理器中 (42)接收其右邻处理器中发来的A的子块 end if (5)if(sqrt(p)}-1)and(mod2=1)and(≠0)then/*其余的奇数号处理器* (5.1)将所存的A的子块在缓冲区buer中做备份 (52)接收其右邻处理器中发来的A的子块 (53)将在缓冲区 buffer中所存的A的子块发送到其左邻处理器中 d if End 当算法执行在√x√P的二维网孔上时,若使用切通CT选路法算法187第2)步的循环 移位时间为2(3+1m"P,第(4)步的单步移位时间为2+1"、运算时间为n3/p。 所以在二维网孔上 Cannon乘法的并行运行时间为Tp=4(s+1yNp+n3/p MPI源程序请参见章末附录 15LU分解 从本小节起我们将对LU分解等矩阵分解的并行计算做一些简单讨论。在许多应用问题 的科学计算中,矩阵的LU分解是基本、常用的一种矩阵运算,它是求解线性方程组的基础, 尤其在解多个同系数阵的线性方程组时特别有用 1.5.1矩阵的LU分解及其串行算法 对于一个n阶非奇异方阵A=c],对A进行LU分解是求一个主对角元素全为1的下三 角方阵L=[与上三角方阵 使A=LU。设A的各阶主子行列式皆非零,U和L的元 素可由下面的递推式求出

算法 18.8 函数 Leftmoveonestep(a)的基本算法 Begin (1)if (j=0) then /*最左端的子块*/ (1.1)将所存的 A 的子块发送到同行最右端子块所在的处理器中 (1.2)接收其右邻处理器中发来的 A 的子块 end if (2)if ((j = sqrt(p)-1) and (j mod 2 = 0)) then /*最右端子块处理器且块列号为偶数*/ (2.1)将所存的 A 的子块发送到其左邻处理器中 (2.2)接收其同行最左端子块所在的处理器发来的 A 的子块 end if (3)if ((j = sqrt(p)-1) and (j mod 2 ≠ 0)) then /*最右端子块处理器且块列号为奇数*/ (3.1)将所存的 A 的子块在缓冲区 buffer 中做备份 (3.2)接收其同行最左端子块所在的处理器发来的 A 的子块 (3.3)将在缓冲区 buffer 中所存的 A 的子块发送到其左邻处理器中 end if (4)if ((j ≠ sqrt(p)-1) and (j mod 2 = 0) and (j ≠ 0)) then /*其余的偶数号处理器*/ (4.1)将所存的 A 的子块发送到其左邻处理器中 (4.2)接收其右邻处理器中发来的 A 的子块 end if (5)if ((j ≠ sqrt(p)-1) and (j mod 2 = 1) and (j ≠ 0)) then /*其余的奇数号处理器*/ (5.1)将所存的 A 的子块在缓冲区 buffer 中做备份 (5.2)接收其右邻处理器中发来的 A 的子块 (5.3)将在缓冲区 buffer 中所存的 A 的子块发送到其左邻处理器中 end if End 当算法执行在 p  p 的二维网孔上时,若使用切通 CT 选路法,算法 18.7 第(2)步的循环 移位时间为 p p n 2(ts tw ) 2 + ,第(4)步的单步移位时间为 p p n 2(ts tw ) 2 + 、运算时间为 n3 / p 。 所以在二维网孔上 Cannon 乘法的并行运行时间为 p n p p n Tp 4(ts tw ) / 3 2 = + + 。 MPI 源程序请参见章末附录。 1.5 LU 分 解 从本小节起我们将对 LU 分解等矩阵分解的并行计算做一些简单讨论。在许多应用问题 的科学计算中，矩阵的 LU 分解是基本、常用的一种矩阵运算，它是求解线性方程组的基础， 尤其在解多个同系数阵的线性方程组时特别有用。 1.5.1 矩阵的 LU 分解及其串行算法 对于一个 n 阶非奇异方阵 A=[aij]，对 A 进行 LU 分解是求一个主对角元素全为 1 的下三 角方阵 L=[lij]与上三角方阵 U=[uij]，使 A=LU。设 A 的各阶主子行列式皆非零，U 和 L 的元 素可由下面的递推式求出：

lik i>k < 在计算过程中,首先计算出U的第一行元素,然后算出L的第一列元素,修改相应A 的元素;再算出U的第二行,L的第二列…,直至算出砌n为止。若一次乘法和加法运算或 次除法运算时间为一个单位时间,则下述LU分解的串行算法189时间复杂度为 (-1)i=(n3-n)/3=O(n) 算法189单处理器上矩阵LU分解串行算法 输入:矩阵 输出:下三角矩阵Lnxn,上三角矩阵Unxn Be (1 ) for k-l to n de (1. 1)for i=k+l to n do a[, k]=a[i, ka[k, k (1.2)for i=k+I to n do for户=k+ I to n do q=q[i小}[*ak nd fo (2 )for i=l to Mi=ali,j1 else u[iF=ali,j (22[i,i]=1 end for 152矩阵LU分解的并行算法 在LU分解的过程中,主要的计算是利用主行i对其余各行j,(>)作初等行变换,各行 计算之间没有数据相关关系,因此可以对矩阵A按行划分来实现并行计算。考虑到在计算 过程中处理器之间的负载均衡,对A采用行交叉划分:设处理器个数为p,矩阵A的阶数为

aij (1)=aij aij (k+1)=aij (k) -likukj         =  = − a u i k i k i k l kk k ik ik ( ) 1 1 0       = a k j k j u k kj kj ( ) 0 在计算过程中，首先计算出 U 的第一行元素，然后算出 L 的第一列元素，修改相应 A 的元素；再算出 U 的第二行，L 的第二列…，直至算出 unn 为止。若一次乘法和加法运算或 一次除法运算时间为一个单位时间，则下述 LU 分解的串行算法 18.9 时间复杂度为 ( 1) ( ) / 3 3 1 i i n n n i  − = − = = O(n 3 )。 算法 18.9 单处理器上矩阵 LU 分解串行算法 输入：矩阵 An×n 输出：下三角矩阵 Ln×n,上三角矩阵 Un×n Begin (1)for k=1 to n do (1.1)for i=k+1 to n do a[i,k]=a[i,k]/a[k,k] end for (1.2)for i=k+1 to n do for j=k+1 to n do a[i,j]=a[i,j]-a[i,k]*a[k,j] end for end for end for (2)for i=1 to n do (2.1)for j=1 to n do if (j<i) then l[i,j]=a[i,j] else u[i,j]=a[i,j] end if end for (2.2)l[i,i]=1 end for End 1.5.2 矩阵 LU 分解的并行算法 在 LU 分解的过程中，主要的计算是利用主行 i 对其余各行 j，(j>i)作初等行变换，各行 计算之间没有数据相关关系，因此可以对矩阵 A 按行划分来实现并行计算。考虑到在计算 过程中处理器之间的负载均衡，对 A 采用行交叉划分：设处理器个数为 p，矩阵 A 的阶数为

n,m=「n/pl,对矩阵A行交叉划分后,编号为=1,…p1)的处理器存有A的第,汁+p… 计+(m-1)行。然后依次以第0.,1,…;n-1行作为主行,将其广播给所有处理器,各处理器利用 主行对其部分行向量做行变换,这实际上是各处理器轮流选出主行并广播。若以编号为 my rank的处理器的第i行元素作为主行,并将它广播给所有处理器,则编号大于等于 my rank的处理器利用主行元素对其第计1,…,m1行数据做行变换,其它处理器利用主行元 素对其第讠…m-1行数据做行变换。具体并行算法框架描述如下 算法18.10矩阵LU分解的并行算法 输入:矩阵Anxn 输出:下三角矩阵Lnxn,上三角矩阵Unxn Begin 对所有处理器 my rank( my rank=0,…,p-l)同时执行如下的算法 for i=0 to m-I do (1)ir( my rank=)then/*当前处理的主行在本处理器* (1.1)叶j/*v为当前处理的主行号* (1.2 )for k-y to n de fk=a{k/*主行元素存到数组∫中* nd fo (1.3)向其它所有处理器广播主行元素 else/*当前处理的主行不在本处理器* (14)=p+ (1.5)接收主行所在处理器广播来的主行元素 (2if( my rank≤then (2. 1 )for k-i+l to m-l do (ia[k, v]=akk, vvI (iife I to n-l do a[k w]=a[k, wl]*a[k,vl end for or (2.2 )for k-i to m-I do (iak, v]=ak, v/fvl (iifor w=v+l to n-l do a[kv=a[kww]°akv or end fo 计算完成后,编号为0的处理器收集各处理器中的计算结果,并从经过初等行变换的矩 阵A中分离出下三角矩阵L和上三角矩阵U。若一次乘法和加法运算或一次除法运算时间 为一个单位时间,则计算时间为(n3-n)/3p;又n-1行数据依次作为主行被广播,通信时间为

n，m = n/ p ，对矩阵 A 行交叉划分后，编号为 i(i=0,1,…,p-1)的处理器存有 A 的第 i, i+p,…, i+(m-1)p 行。然后依次以第 0,1,…,n-1 行作为主行，将其广播给所有处理器，各处理器利用 主行对其部分行向量做行变换，这实际上是各处理器轮流选出主行并广播。若以编号为 my_rank 的处理器的第 i 行元素作为主行，并将它广播给所有处理器，则编号大于等于 my_rank 的处理器利用主行元素对其第 i+1,…,m-1 行数据做行变换，其它处理器利用主行元 素对其第 i,…,m-1 行数据做行变换。具体并行算法框架描述如下： 算法 18.10 矩阵 LU 分解的并行算法 输入：矩阵 An×n 输出：下三角矩阵 Ln×n,上三角矩阵 Un×n Begin 对所有处理器 my_rank(my_rank=0,…,p-1)同时执行如下的算法: for i=0 to m-1 do for j=0 to p-1 do (1)if (my_rank=j) then /*当前处理的主行在本处理器*/ (1.1) v=i*p+j /* v 为当前处理的主行号*/ (1.2)for k=v to n do f[k]=a[i,k] /* 主行元素存到数组 f 中*/ end for (1.3)向其它所有处理器广播主行元素 else /*当前处理的主行不在本处理器*/ (1.4)v=i*p+j (1.5)接收主行所在处理器广播来的主行元素 end if (2)if (my_rank ≤ j) then (2.1)for k=i+1 to m-1 do (i)a[k,v]=a[k,v]/f[v] (ii)for w=v+1 to n-1 do a[k,w]=a[k,w]-f[w]*a[k,v] end for end for else (2.2)for k=i to m-1 do (i)a[k,v]=a[k,v]/f[v]; (ii)for w=v+1 to n-1 do a[k,w]=a[k,w]-f[w]*a[k,v] end for end for end if end for end for End 计算完成后，编号为 0 的处理器收集各处理器中的计算结果，并从经过初等行变换的矩 阵 A 中分离出下三角矩阵 L 和上三角矩阵 U。若一次乘法和加法运算或一次除法运算时间 为一个单位时间，则计算时间为(n 3 -n)/3p；又 n-1 行数据依次作为主行被广播，通信时间为