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

6矩阵运算 矩阵运算是数值计算中最重要的一类运算特别是在线性代数和数值分析中它是一种最 基本的运算。本章讨论的矩阵运算包括矩阵转置、矩阵向量相乘、矩阵乘法、矩阵分解以及 方阵求逆等。在讨论并行矩阵算法时分三步进行:①算法描述及其串行算法:②算法的并行 化及其实现算法框架以及简单的算法分析:③算法实现的MPI源程序,以利于读者实践操 11矩阵转置 111矩阵转置及其串行算法 对于一个n阶方阵A=[o,将其每一下三角元素a(>)沿主对角线与其对称元素a互 换就构成了转置矩阵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 输出:矩阵A1x的转置Anxm 对所有处理器 my rank( my rank=0,…p-1)同时执行如下的算法: (1)计算子块的行号y= 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 交换a[门和叫门 若记t为发送启动时间,l为单位数据传输时间,b为处理器间的延迟时间,则第一步 由于每个子块有πr个元素,又由于通信过程中为了避免死锁,错开下三角子块与上三角 子块的发送顺序,因此子块的交换时间为2(,+1n2/p+1nVP):第二步,假定一对数据 的交换时间为一个单位时间,则局部转置时间为n2/2p。因此所需的并行计算时间 +2,√p+21-+t1√p 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/p1,这些行块依次记为A,A,…,A1,分别存放在标号为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=m2,若处理器个数和向量维数相当,则其时间复杂度为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 0 to n-l do c{=c[i/+a[i]*b[r引 132简单的矩阵并行分块乘法算法 矩阵乘法也可以用分块的思想实现并行,即分块矩阵乘法( Block Matrix Multiplication) 将矩阵A按行划分为P块(P为处理器个数),设u=「m/p,每块含有连续的u行向量,这些 行块依次记为A0A1,…,41,分别存放在标号为0,1,…,p-1的处理器中。对矩阵B按列划分 为P块,记v=[k/p,每块含有连续的v列向量,这些列块依次记为BnB,…Bn1,分别 存放在标号0,1,…p1为的处理器中。将结果矩阵C也相应地同时进行行、列划分,得到 pXP个大小为u×v的子矩阵,记第i行第j列的子矩阵为C,显然有C=A,XB,其中, A1大小为a×n,B大小为n×v 「处理器编号存储内容 存储内容 Ar Bll PI 图12矩阵相乘并行算法中的数据交换

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,…Cφl)。为了避免在通信过程中发生死锁,奇数号及 偶数号处理器的收发顺序被错开,使偶数号处理器先发送后接收;而奇数号处理器先将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 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 计算右邻处理器的标号mpl=( my_ rank+1)modp (4)if(i≠p1)then (4l)ir( my rank mod2=0)then/*编号为偶数的处理器* (i)将所存的B的子块发送到其左邻处理器中 (i)接收其右邻处理器中发来的B的子块 (42if( my rank mod2≠0)then/*编号为奇数的处理器* (i)将所存的B子块在缓冲区 buffer中做备份 (i)接收其右邻处理器中发来的B的子块 (i)将 buffer I所存的B的子块发送到其左邻处理器中 en 设一次乘法和加法运算时间为一个单位时间,由于每个处理器计算p个u×n与n×v 阶的子矩阵相乘,因此计算时间为u考np;所有处理器交换数据p1次,每次的通信量为 γn,通信过程中为了避免死锁,错开奇数号及偶数号处理器的收发顺序,通信时间为 2(p-1)(tx+my*tn)=Onk),所以并行计算时间Tp=p+2p-1)(L+mv)=mnk/p+2(p-1)(+mvtm) 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 源程序请参见所附光盘