6矩阵运算 矩阵运算是数值计算中最重要的一类运算,特别是在线性代数和数值分析中它是一种最 基本的运算。本章讨论的矩阵运算包括矩阵转置、矩阵向量相乘、矩阵乘法、矩阵分解以及 方阵求逆等。在讨论并行矩阵算法时分三步进行:①算法描述及其串行算法;②算法的并行 化及其实现算法框架以及简单的算法分析:③算法实现的MPⅠ源程序,以利于读者实践操 11矩阵转置 111矩阵转置及其串行算法 对于一个n阶方阵A={a],将其每一下三角元素a(>)沿主对角线与其对称元素41互 换就构成了转置矩阵A。假设一对数据的交换时间为一个单位时间,则下述矩阵转置 Matrix Transposing)算法18.1的运行时间为(m2-n)/2=O(m2) 算法18.1单处理器上矩阵转置算法 输入:矩阵Anxn 输出:矩阵Anxn的转置Anxn for i =2 to n do to i-I do 交换叫和a[, endOr endo End 图1.1子块转置 112矩阵转置并行算法 此处主要讨论网孔上的块棋盘划分( Block-Checker Board Partitioning,又称为块状划分) 的矩阵转置算法,它对循环棋盘划分( yclic- Checker Board Partitioning)也同样适用。另外, 超立方上块棋盘划分的矩阵转置算法可参见文献[l 实现矩阵的转置时,若处理器个数为p,且它们的编号依次是0,1,…p1,则将n阶矩
1. 6 矩阵运算 矩阵运算是数值计算中最重要的一类运算,特别是在线性代数和数值分析中,它是一种最 基本的运算。本章讨论的矩阵运算包括矩阵转置、矩阵向量相乘、矩阵乘法、矩阵分解以及 方阵求逆等。在讨论并行矩阵算法时分三步进行：①算法描述及其串行算法；②算法的并行 化及其实现算法框架以及简单的算法分析；③算法实现的 MPI 源程序，以利于读者实践操 作。 1.1 矩阵转置 1.1.1 矩阵转置及其串行算法 对于一个 n 阶方阵 A=[aij]，将其每一下三角元素 aij (i>j)沿主对角线与其对称元素 aji 互 换就构成了转置矩阵 A T。假设一对数据的交换时间为一个单位时间，则下述矩阵转置(Matrix Transposing)算法 18.1 的运行时间为(n 2 -n)/2=O(n 2 )。 算法 18.1 单处理器上矩阵转置算法 输入：矩阵 An×n 输出：矩阵 An×n 的转置 A T n×n Begin for i=2 to n do for j=1 to i-1 do 交换 a[i,j]和 a[j,i] endfor endfor End A41 P12 A42 P13 A43 P14 A12 P1 A13 P2 A14 P3 A11 P0 A21 P4 A23 P6 A24 P7 A22 P5 A31 P8 A32 P9 A34 P11 A33 P10 A44 P15 图 1.1 子块转置 1.1.2 矩阵转置并行算法 此处主要讨论网孔上的块棋盘划分(Block-Checker Board Partitioning，又称为块状划分) 的矩阵转置算法，它对循环棋盘划分(Cyclic-Checker Board Partitioning)也同样适用。另外， 超立方上块棋盘划分的矩阵转置算法可参见文献[1]。 实现矩阵的转置时，若处理器个数为 p，且它们的编号依次是 0,1, …,p-1，则将 n 阶矩
阵A分成p个大小为m×m的子块,m=「n/pp个子块组成一个√x√P的子块阵列。 记其中第i行第j列的子块为Ay,它含有A的第(-1)m+1至第Ⅷm行中的第(-1)m+1至第jm 列的所有元素。对每一处理器按行主方式赋以二维下标,记编号为i的处理器的二维下标为 (),其中=/,=md√,将A的子块存入下标为(a)表示的对应处理器中 这样,转置过程分两步进行:第一步,子块转置,具体过程如图1.1所示:第二步,处理器 内部局部转置。 为了避免对应子块交换数据时处理器发生死锁,可令下三角子块先向与之对应的上三角 子块发送数据,然后从上三角子块接收数据:上三角子块先将数据存放在缓冲区bur中 然后从与之对应的下三角子块接收数据:最后再将缓冲区中的数据发送给下三角子块。具体 并行算法框架描述如下 算法182网孔上的矩阵转置算法 输入:矩阵Axn 输出:矩阵Ax的转置Anx 对所有处理器 my rank( my rank=0…p-1)同时执行如下的算法 (1)计算子块的行号= my ran/ sqrt(p),计算子块的列号l= my rank mod sqrt(p) (2(ux<)then/*对存放下三角块的处理器* (2.1)将所存的子块发送到其对角块所在的处理器中 (2.2)接收其对角块所在的处理器中发来的子块 件对存放上三角块的处理器* (2.3)将所存的子块在缓冲区bur中做备份 (24)接收其对角块所在的处理器中发来的子块 (2.5)将 buffer中所存的子块发送到其对角块所在的处理器中 end if (3)fori=1 to m do/*处理器内部局部转置* forj=l to i do 交换q/和叫 若记t为发送启动时间,l为单位数据传输时间,b为处理器间的延迟时间,则第一步 由于每个子块有πr个元素,又由于通信过程中为了避免死锁,错开下三角子块与上三角 子块的发送顺序,因此子块的交换时间为2(,+1n2/p+1、√P);第二步,假定一对数据 的交换时间为一个单位时间,则局部转置时间为n2/2p。因此所需的并行计算时间 T=+2t,y√p+2 MPI源程序请参见所附光盘
阵 A 分成 p 个大小为 m×m 的子块， m = n / p。p 个子块组成一个 p p 的子块阵列。 记其中第 i 行第 j 列的子块为 Aij，它含有 A 的第(i-1)m+1 至第 im 行中的第(j-1)m+1 至第 jm 列的所有元素。对每一处理器按行主方式赋以二维下标，记编号为 i 的处理器的二维下标为 (v,u)，其中 v = i / p ，u = i mod p ，将 A 的子块存入下标为(v,u)表示的对应处理器中。 这样，转置过程分两步进行：第一步，子块转置，具体过程如图 1.1 所示；第二步，处理器 内部局部转置。 为了避免对应子块交换数据时处理器发生死锁，可令下三角子块先向与之对应的上三角 子块发送数据，然后从上三角子块接收数据；上三角子块先将数据存放在缓冲区 buffer 中， 然后从与之对应的下三角子块接收数据；最后再将缓冲区中的数据发送给下三角子块。具体 并行算法框架描述如下： 算法 18.2 网孔上的矩阵转置算法 输入：矩阵 An×n 输出：矩阵 An×n 的转置 A T n×n Begin 对所有处理器 my_rank(my_rank=0,…,p-1)同时执行如下的算法: (1)计算子块的行号 v=my_rank/ sqrt(p), 计算子块的列号 u=my_rank mod sqrt(p) (2)if (u<v) then /*对存放下三角块的处理器*/ (2.1)将所存的子块发送到其对角块所在的处理器中 (2.2)接收其对角块所在的处理器中发来的子块 else /*对存放上三角块的处理器*/ (2.3)将所存的子块在缓冲区 buffer 中做备份 (2.4)接收其对角块所在的处理器中发来的子块 (2.5)将 buffer 中所存的子块发送到其对角块所在的处理器中 end if (3)for i=1 to m do /*处理器内部局部转置*/ for j=1 to i do 交换 a[i,j]和 a[j,i] end for end for End 若记 ts 为发送启动时间, tw 为单位数据传输时间，th 为处理器间的延迟时间，则第一步 由于每个子块有 n 2 /p 个元素，又由于通信过程中为了避免死锁，错开下三角子块与上三角 子块的发送顺序，因此子块的交换时间为 2( / ) 2 t s + twn p + t h p ；第二步，假定一对数据 的交换时间为一个单位时间，则局部转置时间为 n2 / 2p 。因此所需的并行计算时间 t p p n t p t p n Tp = + s + w + h 2 2 2 2 2 。 MPI 源程序请参见所附光盘
12矩阵向量乘法 121矩阵向量乘法及其串行算法 矩阵向量乘法 Matrix-Vector Multiplication)是将一个n×n阶方阵A=a]乘以n×1的向 量B={b,b2…,b]得到一个具有n个元素的列向量C-C,c,…c假设一次乘法和加法 运算时间为一个单位时间,则下述矩阵向量乘法算法83的时间复杂度为Omn3) 算法183单处理器上矩阵向量乘算法 输入:Anxn,Bnx 输出:Cn×1 for i=0 to n-l do c[=0 forj=0 to n-l do cFc+aibl end for end for 122矩阵向量乘法的并行算法 矩阵向量乘法同样可以有带状划分( Striped Partitioning)和棋盘划分( Checker board Partitioning)两种并行算法。以下仅讨论行带状划分矩阵向量乘法,列带状划分矩阵向量乘法 是类似的。设处理器个数为p,对矩阵A按行划分为p块,每块含有连续的m行向量 m=「n/p,这些行块依次记为A,A,…,An,分别存放在标号为0,1…p1的处理器中, 同时将向量B广播给所有处理器。各处理器并行地对存于局部数组a中的行块A和向量B 做乘积操作,具体并行算法框架描述如下: 算法184行带状划分的矩阵向量乘并行算法 输入:Anxn,Bnx1 输出:C nxI Be 对所有处理器 my rank( my rank=0,…,p-1)同时执行如下的算法 for i=0 to m-1 do c()=0.0 for户=0ton-ldo c{=c可+a[]*b end for end fe End 假设一次乘法和加法运算时间为一个单位时间,不难得出行带状划分的矩阵-向量乘算 法184并行计算时间T=m/,若处理器个数和向量维数相当,则其时间复杂度为O(n) MPI源程序请参见所附光盘
1.2 矩阵-向量乘法 1.2.1 矩阵-向量乘法及其串行算法 矩阵-向量乘法(Matrix-Vector Multiplication)是将一个 n×n 阶方阵 A=[aij]乘以 n×1 的向 量 B=[b1,b2, …,bn] T得到一个具有 n 个元素的列向量 C=[c1,c2, …,cn] T。假设一次乘法和加法 运算时间为一个单位时间，则下述矩阵向量乘法算法 18.3 的时间复杂度为 O(n 2 )。 算法 18.3 单处理器上矩阵-向量乘算法 输入：An×n，Bn×1 输出：Cn×1 Begin for i=0 to n-1 do c[i]= 0 for j=0 to n-1 do c[i]=c[i] + a[i,j]*b[j] end for end for End 1.2.2 矩阵-向量乘法的并行算法 矩阵-向量乘法同样可以有带状划分(Striped Partitioning)和棋盘划分(Checker Board Partitioning)两种并行算法。以下仅讨论行带状划分矩阵-向量乘法,列带状划分矩阵-向量乘法 是类似的。设处理器个数为 p，对矩阵 A 按行划分为 p 块，每块含有连续的 m 行向量， m = n / p ，这些行块依次记为 A0, A1, …, Ap-1，分别存放在标号为 0,1,…,p-1 的处理器中， 同时将向量 B 广播给所有处理器。各处理器并行地对存于局部数组 a 中的行块 Ai 和向量 B 做乘积操作，具体并行算法框架描述如下： 算法 18.4 行带状划分的矩阵-向量乘并行算法 输入：An×n，Bn×1 输出：Cn×1 Begin 对所有处理器 my_rank(my_rank=0,…, p-1)同时执行如下的算法: for i=0 to m-1 do c(i)=0.0 for j=0 to n-1 do c[i] = c[i]+a[i,j]*b[j] end for end for End 假设一次乘法和加法运算时间为一个单位时间，不难得出行带状划分的矩阵-向量乘算 法 18.4 并行计算时间 Tp=n 2 /p，若处理器个数和向量维数相当，则其时间复杂度为 O(n)。 MPI 源程序请参见所附光盘
13行列划分矩阵乘法 个m×n阶矩阵A=[a乘以一个n×k的矩阵B=[b就可以得到一个m×k的矩阵 C=[c],它的元素c为A的第i行向量与B的第j列向量的内积。矩阵相乘的关键是相乘的 两个元素的下标要满足一定的要求(即对准)。为此常采用适当旋转矩阵元素的方法(如后面将 要阐述的 Cannon乘法),或采用适当复制矩阵元素的办法(如DNS算法),或采用流水线的 办法使元素下标对准,后两种方法读者可参见文献[]l 131矩阵相乘及其串行算法 由矩阵乘法定义容易给出其串行算法18.5,若一次乘法和加法运算时间为一个单位时 间,则显然其时间复杂度为O(mk) 算法18.5单处理器上矩阵相乘算法 输入:Amxn,Bn×k 输出:C for i=0 to m-1 do for j=0 to k-l do for r0 to n- do c{=c[i/+a[i]*b[r引 end for 132简单的矩阵并行分块乘法算法 矩阵乘法也可以用分块的思想实现并行,即分块矩阵乘法( Block Matrix Multiplication) 将矩阵A按行划分为P块(p为处理器个数),设u=「m/p,每块含有连续的u行向量,这些 行块依次记为A0A1,…41,分别存放在标号为0,1,…,p-1的处理器中。对矩阵B按列划分 为P块,记v=[k/p,每块含有连续的v列向量,这些列块依次记为BB1,…Bn1,分别 存放在标号0,1,…p1为的处理器中。将结果矩阵C也相应地同时进行行、列划分,得到 pXP个大小为u×v的子矩阵,记第i行第j列的子矩阵为C,显然有C=A,XB,其中, A1大小为a×n,B大小为n×v 「处理器编号存储内容 存储内容 Ar Bll PI 图1.2矩阵相乘并行算法中的数据交换
1.3 行列划分矩阵乘法 一个 m×n 阶矩阵 A=[aij]乘以一个 n×k 的矩阵 B=[bij]就可以得到一个 m×k 的矩阵 C=[cij]，它的元素 cij 为 A 的第 i 行向量与 B 的第 j 列向量的内积。矩阵相乘的关键是相乘的 两个元素的下标要满足一定的要求(即对准)。为此常采用适当旋转矩阵元素的方法(如后面将 要阐述的 Cannon 乘法)，或采用适当复制矩阵元素的办法(如 DNS 算法)，或采用流水线的 办法使元素下标对准，后两种方法读者可参见文献[1]。 1.3.1 矩阵相乘及其串行算法 由矩阵乘法定义容易给出其串行算法 18.5，若一次乘法和加法运算时间为一个单位时 间，则显然其时间复杂度为 O(mnk) 。 算法 18.5 单处理器上矩阵相乘算法 输入：Am×n，Bn×k 输出：Cm×k Begin for i=0 to m-1 do for j=0 to k-1 do c[i,j]=0 for r=0 to n-1 do c[i,j]= c[i,j]+a[i,r]*b[r,j] end for end for end for End 1.3.2 简单的矩阵并行分块乘法算法 矩阵乘法也可以用分块的思想实现并行，即分块矩阵乘法(Block Matrix Multiplication)， 将矩阵 A 按行划分为 p 块(p 为处理器个数)，设 u = m/ p ，每块含有连续的 u 行向量，这些 行块依次记为 A0,A1, …,Ap-1，分别存放在标号为 0,1,…, p-1 的处理器中。对矩阵 B 按列划分 为 P 块，记 v = k / p ，每块含有连续的 v 列向量，这些列块依次记为 B0,B1, …,Bp-1，分别 存放在标号 0,1, …, p-1 为的处理器中。将结果矩阵 C 也相应地同时进行行、列划分，得到 p×p 个大小为 u×v 的子矩阵，记第 i 行第 j 列的子矩阵为 Cij，显然有 Cij= Ai×Bj，其中， A i 大小为 u×n，Bj 大小为 n×v。 A0 处理器编号 存储内容 B0 0 A1 B1 A2 B2 Ap-1 Bp-1 1 2 P-1 … … … (a) A0 处理器编号 存储内容 B1 0 A1 B2 A2 B3 Ap-1 B0 1 2 P-1 … … … (b) 图 1.2 矩阵相乘并行算法中的数据交换
开始,各处理器的存储内容如图12(a)所示。此时各处理器并行计算C=AXB其中 产=0,1,…p-1,此后第i号处理器将其所存储的B的列块送至第μ-1号处理器(第0号处理器 将B的列块送至第p-1号处理器中,形成循环传送),各处理器中的存储内容如图1.2(b)所 示。它们再次并行计算C=A1XB,这里产=(+1)modp。B的列块在各处理器中以这样的方 式循环传送p-1次并做p次子矩阵相乘运算,就生成了矩阵C的所有子矩阵。编号为i的处 理器的内部存储器存有子矩阵Co,C,…Cl)。为了避免在通信过程中发生死锁,奇数号及 偶数号处理器的收发顺序被错开,使偶数号处理器先发送后接收;而奇数号处理器先将B 的列块存于缓冲区bur中,然后接收编号在其后面的处理器所发送的B的列块,最后再将 缓冲区中原矩阵B的列块发送给编号在其前面的处理器,具体并行算法框架描述如下: 算法186矩阵并行分块乘法算法 输入:Amxn,Bnx 输出:Cm 对所有处理器 my rank( my rank=0,…p-1)同时执行如下的算法: (1)目前计算C的子块号l=(+ my rank)modp (2)for =0 to u-I do for产0 to l-I do c{.==0 for so to n-l do c{.二小=c[,]+a[=,s]*b|s end for end for (3)计算左邻处理器的标号mml=(p+ my rank-1)modp 计算右邻处理器的标号mp1=( my_ rank+1)modp (4)if(i≠p1)then (4l)ir( my rank mod2=0)then/*编号为偶数的处理器*/ (i)将所存的B的子块发送到其左邻处理器中 (i)接收其右邻处理器中发来的B的子块 end if (42if( my rank mod2≠0)then/*编号为奇数的处理器* (i)将所存的B子块在缓冲区buer中做备份 (i)接收其右邻处理器中发来的B的子块 (i)将 buffer I所存的B的子块发送到其左邻处理器中 en 设一次乘法和加法运算时间为一个单位时间,由于每个处理器计算p个u×n与n×v 阶的子矩阵相乘,因此计算时间为u市p;所有处理器交换数据p1次,每次的通信量为 γn,通信过程中为了避免死锁,错开奇数号及偶数号处理器的收发顺序,通信时间为 2(P1)(1+m*m)=O(nk),所以并行计算时间Tp=mp+2p1)+mtm)=mnk/p+2(p1)(4+mtm) MPI源程序请参见所附光盘
开始，各处理器的存储内容如图 1.2 (a)所示。此时各处理器并行计算 Cii= Ai×Bj 其中 i=0,1,…,p-1，此后第 i 号处理器将其所存储的 B 的列块送至第 i-1 号处理器（第 0 号处理器 将 B 的列块送至第 p-1 号处理器中，形成循环传送），各处理器中的存储内容如图 1.2 (b)所 示。它们再次并行计算 Cij= A i×B j，这里 j=(i+1)modp。B 的列块在各处理器中以这样的方 式循环传送 p-1 次并做 p 次子矩阵相乘运算，就生成了矩阵 C 的所有子矩阵。编号为 i 的处 理器的内部存储器存有子矩阵 Ci0,Ci1,…,Ci(p-1)。为了避免在通信过程中发生死锁，奇数号及 偶数号处理器的收发顺序被错开，使偶数号处理器先发送后接收；而奇数号处理器先将 B 的列块存于缓冲区 buffer 中，然后接收编号在其后面的处理器所发送的 B 的列块，最后再将 缓冲区中原矩阵 B 的列块发送给编号在其前面的处理器，具体并行算法框架描述如下： 算法 18.6 矩阵并行分块乘法算法 输入：Am×n, Bn×k, 输出：Cm×k Begin 对所有处理器 my_rank(my_rank=0,…,p-1)同时执行如下的算法: (1)目前计算 C 的子块号 l=(i+my_rank) mod p (2)for z=0 to u-1 do for j=0 to v-1 do c[l,z,j]=0 for s=0 to n-1 do c[l,z,j]=c[l,z,j]+ a[z,s]*b[s,j] end for end for end for (3)计算左邻处理器的标号 mm1=(p+my_rank-1) mod p 计算右邻处理器的标号 mp1=(my_rank+1) mod p (4)if (i≠p-1) then (4.1)if (my_rank mod 2= 0) then /*编号为偶数的处理器*/ (i)将所存的 B 的子块发送到其左邻处理器中 (ii)接收其右邻处理器中发来的 B 的子块 end if (4.2)if (my_rank mod 2≠ 0) then /*编号为奇数的处理器*/ (i)将所存的 B 子块在缓冲区 buffer 中做备份 (ii)接收其右邻处理器中发来的 B 的子块 (iii)将 buffer 中所存的 B 的子块发送到其左邻处理器中 end if end if End 设一次乘法和加法运算时间为一个单位时间，由于每个处理器计算 p 个 u×n 与 n×v 阶的子矩阵相乘，因此计算时间为 u*v*n*p；所有处理器交换数据 p-1 次，每次的通信量为 v*n，通信过程中为了避免死锁，错开奇数号及偶数号处理器的收发顺序，通信时间为 2(p-1)(ts+nv*tw)=O(nk)，所以并行计算时间 Tp=uvnp+2(p-1)(ts+nvtw)=mnk / p+2(p-1)(ts+nvtw)。 MPI 源程序请参见所附光盘
14 Cannon乘法 141 Cannon乘法的原理 Cannon算法是一种存储有效的算法。为了使两矩阵下标满足相乘的要求,它和上一节 的并行分块乘法不同,不是仅仅让B矩阵的各列块循环移动,而是有目的地让A的各行块 以及B的各列块皆施行循环移位,从而实现对C的子块的计算。将矩阵A和B分成p个方 块A和Bn,(0≤,j≤VP-1,每块大小为「mP1x「m√F,并将它们分配给√x√P个处 理器(Po,P1…,P--1)开始时处理器P存放块Ay和Bn,并负责计算块Cn,然 后算法开始执行: (1)将块A向左循环移动i步;将块B向上循环移动j步; (2)P执行乘加运算后将块A向左循环移动1步,块Bn向上循环移动1步 (3)重复第2)步,总共执行√P次乘加运算和√P次块A和Bn的循环单步移位。 142 Cannon乘法的并行算法 图13示例了在16个处理器上,用 Cannon算法执行A×x4×B4x4的过程。其中(a)和(b) 对应于上述算法的第(1)步;(c)、(d)、(e)、(f)对应于上述算法的第(2)和第(3)步。在算法第() 步时,A矩阵的第0列不移位,第Ⅰ行循环左移1位,第2行循环左移2位,第3行循环左 移3位:类似地,B矩阵的第0行不移位,第1列循环上移Ⅰ位,第2列循环上移2列,第 3列循环上移3列。这样 Cannon算法具体描述如下: 算法187 Cannon乘法算法 输入:Anxn,Bnxn 输出:Cnxn 对所有处理器 my_ rank( my rank=0,…p-1)同时执行如下的算法 (1)计算子块的行号r= my rank/sqrt(p) 计算子块的列号产 my rank mod sqrt(p) (2)fork=0to√p-ldo if(pk) then Leftmoveonestep(a) end if陳*a循环左移至同行相邻处理器中* if(>k) then Upmoveonestep(b) end if*b循环上移至同列相邻处理器中*/ (3for i=0 to m-l do for户=0tom-ldo c[/]=0 end for end for (4 )for k=0 to p-l de for i=0 to m-1 do or=0 to m-l do
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[i/+a[,k1]*b[k1 Leftmoveonestep(a)/子块a循环左移至同行相邻的处理器中* Upmoveonestep(b)/*子块b循环上移至同列相邻的处理器中*/ nd fo End A3 (a)阵起始对准 b)B阵起始对准 (c)对准后的A和B (d)第一次移位后的子阵位置 B3 303 ::2:2:签2:2 (e)第二次移位后的子阵位置 (f)第三次移位后的子阵位置 图1.316个处理器上的 Cannon乘法过程 这里函数 Leftmoveonestep(a)表示子块a在编号处于同行的处理器之间以循环左移的方 式移至相邻的处理器中;函数 Upmoveonestepl(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中做备份 (32)接收其同行最左端子块所在的处理器发来的A的子块 (33)将在缓冲区 buffer中所存的A的子块发送到其左邻处理器中 end if (4)if(≠sqpt(p)-)and(md2=0)and(≠0)then/其余的偶数号处理器 (4.1)将所存的A的子块发送到其左邻处理器中 (42)接收其右邻处理器中发来的A的子块 end if (5)if(sqrt(p)-1)and(mod2=1)and(≠0)then/*其余的奇数号处理器* (5.1)将所存的A的子块在缓冲区 buffer中做备份 (52)接收其右邻处理器中发来的A的子块 (53)将在缓冲区 buffer中所存的A的子块发送到其左邻处理器中 end if End 当算法执行在√x√P的二维网孔上时若使用切通CT选路法算法187第2)步的循环 移位时间为2(3+1m"-P,第(4)步的单步移位时间为2+1"、运算时间为n3/p。 所以在二维网孔上 Cannon乘法的并行运行时间为Tp=4(s+1m")p+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 )作初等行变换,各行 计算之间没有数据相关关系,因此可以对矩阵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 (ji)作初等行变换，各行 计算之间没有数据相关关系，因此可以对矩阵 A 按行划分来实现并行计算。考虑到在计算 过程中处理器之间的负载均衡，对 A 采用行交叉划分：设处理器个数为 p，矩阵 A 的阶数为
n,m=「n/p,对矩阵A行交叉划分后,编号为(=0.1,…p1)的处理器存有A的第,计p…, 计+(m-1)行。然后依次以第0.,1,…;n-1行作为主行,将其广播给所有处理器,各处理器利用 主行对其部分行向量做行变换,这实际上是各处理器轮流选出主行并广播。若以编号为 my rank的处理器的第i行元素作为主行,并将它广播给所有处理器,则编号大于等于 my rank的处理器利用主行元素对其第汁1;…m-1行数据做行变换,其它处理器利用主行元 素对其第讠…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]=a[k, vl/vI (ii)fe I to n-l do a[k w]=a[k, w]fw]*alk end for or (2.2 )for k-i to m-I do (ia[k, v]=a[k, vlv a[kv=a[kww]°akv or 计算完成后,编号为0的处理器收集各处理器中的计算结果,并从经过初等行变换的矩 阵A中分离出下三角矩阵L和上三角矩阵U。若一次乘法和加法运算或一次除法运算时间 为一个单位时间,则计算时间为(m-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 行数据依次作为主行被广播，通信时间为