6 矩阵运算
矩阵运算是数值计算中最重要的一类运算,特别是在线性代数和数值分析中,它是一种最基本的运算。本章讨论的矩阵运算包括矩阵转置、矩阵向量相乘、矩阵乘法、矩阵分解以及方阵求逆等。在讨论并行矩阵算法时分三步进行：①算法描述及其串行算法；②算法的并行化及其实现算法框架以及简单的算法分析；③算法实现的MPI源程序，以利于读者实践操作。
矩阵转置
矩阵转置及其串行算法
对于一个n阶方阵A=[aij]，将其每一下三角元素aij (i>j)沿主对角线与其对称元素aji互换就构成了转置矩阵AT。假设一对数据的交换时间为一个单位时间，则下述矩阵转置(Matrix Transposing)算法18.1的运行时间为(n2-n)/2=O(n2)。
算法18.1 单处理器上矩阵转置算法
输入：矩阵An×n
输出：矩阵An×n的转置ATn×n
Begin
for i=2 to n do
for j=1 to i-1 do
交换a[i,j]和a[j,i]
endfor
endfor
End
图 18.1 子块转置
矩阵转置并行算法
此处主要讨论网孔上的块棋盘划分(Block-Checker Board Partitioning，又称为块状划分)的矩阵转置算法，它对循环棋盘划分(Cyclic-Checker Board Partitioning)也同样适用。另外，超立方上块棋盘划分的矩阵转置算法可参见文献[1]。
实现矩阵的转置时，若处理器个数为p，且它们的编号依次是0,1, …,p-1，则将n阶矩阵A分成p个大小为m×m的子块，。p个子块组成一个的子块阵列。记其中第i行第j列的子块为Aij，它含有A的第(i-1)m+1至第im行中的第(j-1)m+1至第jm列的所有元素。对每一处理器按行主方式赋以二维下标，记编号为i的处理器的二维下标为(v,u)，其中，，将A的子块存入下标为(v,u)表示的对应处理器中。
这样，转置过程分两步进行：第一步，子块转置，具体过程如图 18.1所示；第二步，处理器内部局部转置。
为了避免对应子块交换数据时处理器发生死锁，可令下三角子块先向与之对应的上三角子块发送数据，然后从上三角子块接收数据；上三角子块先将数据存放在缓冲区buffer中，然后从与之对应的下三角子块接收数据；最后再将缓冲区中的数据发送给下三角子块。具体并行算法框架描述如下：
算法18.2 网孔上的矩阵转置算法
输入：矩阵An×n
输出：矩阵An×n的转置ATn×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为处理器间的延迟时间，则第一步由于每个子块有n2/p个元素，又由于通信过程中为了避免死锁，错开下三角子块与上三角子块的发送顺序，因此子块的交换时间为；第二步，假定一对数据的交换时间为一个单位时间，则局部转置时间为。因此所需的并行计算时间。
MPI源程序请参见所附光盘。
矩阵-向量乘法
矩阵-向量乘法及其串行算法
矩阵-向量乘法(Matrix-Vector Multiplication)是将一个n×n阶方阵A=[aij]乘以n×1的向量B=[b1,b2, …,bn]T得到一个具有n个元素的列向量C=[c1,c2, …,cn]T。假设一次乘法和加法运算时间为一个单位时间，则下述矩阵向量乘法算法18.3的时间复杂度为O(n2)。
算法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
矩阵-向量乘法的并行算法
矩阵-向量乘法同样可以有带状划分(Striped Partitioning)和棋盘划分(Checker Board Partitioning)两种并行算法。以下仅讨论行带状划分矩阵-向量乘法,列带状划分矩阵-向量乘法是类似的。设处理器个数为p，对矩阵A按行划分为p块，每块含有连续的m行向量，，这些行块依次记为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=n2/p，若处理器个数和向量维数相当，则其时间复杂度为O(n)。
MPI源程序请参见所附光盘。
行列划分矩阵乘法
一个m×n阶矩阵A=[aij]乘以一个n×k的矩阵B=[bij]就可以得到一个m×k的矩阵C=[cij]，它的元素cij为A的第i行向量与B的第j列向量的内积。矩阵相乘的关键是相乘的两个元素的下标要满足一定的要求(即对准)。为此常采用适当旋转矩阵元素的方法(如后面将要阐述的Cannon乘法)，或采用适当复制矩阵元素的办法(如DNS算法)，或采用流水线的办法使元素下标对准，后两种方法读者可参见文献[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
简单的矩阵并行分块乘法算法
矩阵乘法也可以用分块的思想实现并行，即分块矩阵乘法(Block Matrix Multiplication)，将矩阵A按行划分为p块(p为处理器个数)，设，每块含有连续的u行向量，这些行块依次记为A0,A1, …,Ap-1，分别存放在标号为0,1,…, p-1的处理器中。对矩阵B按列划分为P块，记，每块含有连续的v列向量，这些列块依次记为B0,B1, …,Bp-1，分别存放在标号0,1, …, p-1为的处理器中。将结果矩阵C也相应地同时进行行、列划分，得到p×p个大小为u×v的子矩阵，记第i行第j列的子矩阵为Cij，显然有Cij= A i×B j，其中，A i 大小为u×n，B j大小为n×v。
图 18.2 矩阵相乘并行算法中的数据交换
开始，各处理器的存储内容如图 18.2 (a)所示。此时各处理器并行计算Cii= A i×B j其中i=0,1,…,p-1，此后第i号处理器将其所存储的B的列块送至第i-1号处理器（第0号处理器将B的列块送至第p-1号处理器中，形成循环传送），各处理器中的存储内容如图 18.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源程序请参见所附光盘。
Cannon乘法
Cannon乘法的原理
Cannon算法是一种存储有效的算法。为了使两矩阵下标满足相乘的要求，它和上一节的并行分块乘法不同，不是仅仅让B矩阵的各列块循环移动，而是有目的地让A的各行块以及B的各列块皆施行循环移位，从而实现对C的子块的计算。将矩阵A和B分成p个方块Aij和Bij，，每块大小为，并将它们分配给个处理器。开始时处理器P ij存放块Aij和Bij，并负责计算块Cij，然后算法开始执行：
⑴将块Aij向左循环移动i步；将块Bij向上循环移动j步；
⑵Pij执行乘加运算后将块Aij向左循环移动1步，块Bij向上循环移动1步；
⑶重复第⑵步，总共执行次乘加运算和次块Aij和Bij的循环单步移位。
Cannon乘法的并行算法
图 18.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 -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 -1 do
for i=0 to m-1 do
for j=0 to m-1 do
for k1=0 to m-1 do
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
图 18.3 16个处理器上的Cannon乘法过程
这里函数Leftmoveonestep(a)表示子块a在编号处于同行的处理器之间以循环左移的方式移至相邻的处理器中；函数Upmoveonestep(b)表示子块b在编号处于同列的处理器之间以循环上移的方式移至相邻的处理器中。这里我们以函数Leftmoveonestep(a)为例，给出处理器间交换数据的过程：
算法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
当算法执行在的二维网孔上时,若使用切通CT选路法,算法18.7第(2)步的循环移位时间为,第(4)步的单步移位时间为、运算时间为。所以在二维网孔上Cannon乘法的并行运行时间为。
MPI源程序请参见章末附录。
LU 分 解
从本小节起我们将对LU分解等矩阵分解的并行计算做一些简单讨论。在许多应用问题的科学计算中，矩阵的LU分解是基本、常用的一种矩阵运算，它是求解线性方程组的基础，尤其在解多个同系数阵的线性方程组时特别有用。
矩阵的LU分解及其串行算法
对于一个n阶非奇异方阵A=[aij]，对A进行LU分解是求一个主对角元素全为1的下三角方阵L=[lij]与上三角方阵U=[uij]，使A=LU。设A的各阶主子行列式皆非零，U和L的元素可由下面的递推式求出：
aij(1)=aij
aij(k+1)=aij(k)-likukj
在计算过程中，首先计算出U的第一行元素，然后算出L的第一列元素，修改相应A的元素；再算出U的第二行，L的第二列…，直至算出unn为止。若一次乘法和加法运算或一次除法运算时间为一个单位时间，则下述LU分解的串行算法18.9时间复杂度为= O(n3)。
算法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
矩阵LU分解的并行算法
在LU分解的过程中，主要的计算是利用主行i对其余各行j，(j>i)作初等行变换，各行计算之间没有数据相关关系，因此可以对矩阵A按行划分来实现并行计算。考虑到在计算过程中处理器之间的负载均衡，对A采用行交叉划分：设处理器个数为p，矩阵A的阶数为n，，对矩阵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。若一次乘法和加法运算或一次除法运算时间为一个单位时间，则计算时间为(n3-n)/3p；又n-1行数据依次作为主行被广播，通信时间为(n-1)(ts+ntw)logp=O(n log p)，因此并行计算时间Tp=(n3-n)/3p+(n-1)(ts+ntw)logp。
MPI源程序请参见章末附录。
QR 分 解
A=[aij] 为一个n阶实矩阵，对A进行QR分解，就是求一个非奇异(Nonsingular)方阵Q与上三角方阵R，使得A=QR。其中方阵Q满足：QT=Q?1，称为正交矩阵(Orthogonal Matrix)，因此QR分解又称为正交三角分解。
矩阵QR分解的串行算法
对A进行正交三角分解，可用一系列平面旋转矩阵左乘A，以使A的下三角元素逐个地被消为0。设要消去的下三角元素为aij(i>j)，则所用的平面旋转方阵为：
其中，除了(i,i)(j,j)元素等于c，(i,j)元素与(j,i)元素分别等于?s与s以外，其余元素皆与单位方阵I的对应元素相等，而c,s由下式计算：
其中，θ为旋转角度。这样，在旋转后所得到的方阵中，元素为0，与A相比仅在第i行、第行上的元素不同：
消去A的r=n(n-1)/2个下三角元素后得到上三角阵R，实际上是用r个旋转方阵Q1,Q2,…,Qr左乘A，即：
R= Qr Qr-1…Q1A
设Q0=Qr Qr-1…Q1，易知Q0为一正交方阵，则有：
R= Q0A
即 A = Q0?1R= Q0TR= QR
其中Q= Q0?1 =Q0T为一正交方阵，而Q0可以通过对单位阵I进行同样的变换而得到, 这样可得到A的正交三角分解。QR分解串行算法如下：
算法18.11 单处理器上矩阵的QR分解串行算法算法
输入：矩阵An×n，单位矩阵Q
输出：矩阵Qn×n，矩阵 Rn×n
Begin
(1)for j=1 to n do
for i=j+1 to n do
(i)sq=sqrt(a[j,j]*a[j,j]+a[i,j]*a[i,j])
c= a[j,j]/sq
s= a[i,j]/sq
(ii)for k=1 to n do
aj[k]= c*a[j,k] + s*a[i,k]
qj[k]=c*q[j,k] + s*q[i,k]
ai[k]= ?s*a[j,k] + c*a[i,k]
qi[k]= ?s*q[j,k] + c*q[i,k]
end for
(iii)for k=1 to n do
a[j,k]=aj[k]
q[j,k]=qj[k]
a[i,k]=ai[k]
q[i,k]=qi[k]
end for
end for
end for
(2)R=A
(3)MATRIX_TRANSPOSITION(Q) /* 对Q实施算法18.1的矩阵转置*/
End
算法18.11要对n(n-1)/2个下三角元素进行消去，每消去一个元素要对两行元素进行旋转变换，而对一个元素进行旋转变换又需要4次乘法和2次加法，所以若取一次乘法或一次加法运算时间为一个单位时间，则该算法的计算时间为n(n-1)/2*2n*6=6n2(n-1)=O(n3)。
矩阵QR分解的并行算法
由于QR分解中消去aij时，同时要改变第i行及第j行两行的元素，而在LU分解中，仅利用主行i(i<j)变更第j行的元素。因此QR分解并行计算中对数据的划分与分布与LU分解就不一样。设处理器个数为p，对矩阵A按行划分为p块，每块含有连续的m行向量，，这些行块依次记为A0,A1, …,Ap-1，分别存放在标号为0,1,…,p-1的处理器中。
在0号处理器中，计算开始时，以第0行数据作为主行，依次和第0,1,…,m-1行数据做旋转变换，计算完毕将第0行数据发送给1号处理器，以使1号处理器将收到的第0行数据作为主行和自己的m行数据依次做旋转变换；在此同时，0号处理器进一步以第1行数据作为主行，依次和第2,3,…,m-1行数据做旋转变换，计算完成后将第1行数据发送给1号处理器；如此循环下去。一直到以第m-1行数据作为主行参与旋转变换为止。
在1号处理器中，首先以收到的0号处理器的第0行数据作为主行和自己的m行数据做旋转变换，计算完将该主行数据发送给2号处理器；然后以收到的0号处理器的第1行数据作为主行和自己的m行数据做旋转变换，再将该主行数据发送给2号处理器；如此循环下去，一直到以收到的0号处理器的第m-1行数据作为主行和自己的m行数据做旋转变换并将第m-1行数据发送给2号处理器为止。然后，1号处理器以自己的第0行数据作为主行，依次和第1, 2,…,m-1行数据做旋转变换，计算完毕将第0行数据发送给2号处理器；接着以第1行数据作为主行，依次和第2,3,…,m-1行数据做旋转变换，计算完毕将第1行数据发送给2号处理器；如此循环下去。一直到以第m-1行数据作为主行参与旋转变换为止。除了p-1号处理器以外的所有处理器都按此规律先顺序接收前一个处理器发来的数据，并作为主行和自己的m行数据作旋转变换，计算完毕将该主行数据发送给后一个处理器。然后依次以自己的m行数据作主行对主行后面的数据做旋转变换，计算完毕将主行数据发给后一个处理器。
在p-1号处理器中，先以接收到的数据作为主行参与旋转变换，然后依次以自己的m行数据作为主行参与旋转变换。所不同的是，p-1号处理器作为最后一个处理器，负责对经过计算的全部数据做汇总。具体并行算法框架描述如下：
算法18.12 矩阵QR分解并行算法
输入：矩阵An×n，单位矩阵Q
输出：矩阵Qn×n，矩阵 Rn×n
Begin
对所有处理器my_rank(my_rank=0,…, p-1)同时执行如下的算法:
(1)if (my_rank=0) then /*0号处理器*/
(1.1)for j=0 to m-2 do
(i)for i=j+1 to m-1 do
Turnningtransform( ) /*旋转变换*/
end for
(ii)将旋转变换后的A和Q的第j 行传送到第1号处理器
end for
(1.2)将旋转变换后的A和Q的第m-1行传送到第1号处理器
end if
(2)if ((my_rank>0) and (my_rank<(p-1))) then /*中间处理器*/
(2.1) for j=0 to my_rank*m-1 do
(i)接收左邻处理器传送来的A和Q子块中的第j 行
(ii)for i=0 to m-1 do
Turnningtransform( ) /*旋转变换*/
end for
(iii)将旋转变换后的A和Q子块中的第j 行传送到右邻处理器
end for
(2.2) for j=0 to m-2 do
(i)z=my_rank*m
(ii)for i=j+1 to m-1 do
Turnningtransform( ) /*旋转变换*/
end for
(iii)将旋转变换后的A和Q子块中的第j 行传送到右邻处理器
end for
(2.3)将旋转变换后的A和Q子块中的第m-1行传送到右邻处理器
end if
(3)if (my_rank= (p-1)) then /*p-1号处理器*/
(3.1) for j=0 to my_rank*m-1 do
(i)接收左邻处理器传送来的A和Q子块中的第j 行
(ii)for i=0 to m-1 do
Turnningtransform( ) /*旋转变换*/
end for
end for
(3.2)for j=0 to m-2 do
for i=j+1 to m-1 do f
Turnningtransform( ) /*旋转变换*/
end for
end for
end if
End
当所有的计算完成后，编号为p-1的处理器负责收集各处理器中的计算结果，R矩阵为经旋转变换的A矩阵，Q矩阵为经旋转变换的Q的转置。QR分解并行计算的时间复杂度分析因每个处理器的开始计算时间不同而不同于其它算法，每个处理器都要顺序对其局部存储器中的m行向量做自身的旋转变换，其时间用Tcomputer表示。以16阶矩阵为例，4个处理器对其进行QR分解的时间分布图如图 18.4所示，这里i(j)表示以接收到的第j号处理器的第i行数据作为主行和自己的m行数据做旋转变换。?t表示第p-1号处理器与第0号处理器开始计算时间的间隔。
图 18.4 QR分解并行计算的时间分布图
此处不妨以第p-1号处理器为对象分析并行计算时间，若取一次乘法或一次加法运算时间为一个单位时间，由于对一个元素进行旋转变换需要4次乘法和2次加法，因而需要六个单位时间。要对m(m-1)/2个下三角元素进行消去，每消去一个元素要对两行元素进行旋转变换，故自身的旋转变换时间Tcomputer =6mn(m-1)；第p-1号处理器依次接收前n-m行数据作为主行参与旋转变换，其计算时间为(n-m)*m*2n*6=12nm(n-m)，通信时间为2(n-m)(ts+ ntw)；时间间隔?t即第p-1号处理器等待第0号处理器第0行数据到达的时间，第0行数据在第0号处理器中参与(m-1)行旋转变换，然后被发送给第1号处理器，这段计算时间为(m-1)*2n*6，通信时间为2(ts+ ntw)；接着第0行数据经过中间(p-2)个处理器的旋转变换被发送给第p-1号处理器，这段计算时间为(p-2)m*2n*6，通信时间为4(p-2)(ts+ntw)，所以?t=12n(mp-m-1)+(4p-6)(ts+ntw), 并行计算时间为Tp=12n2-18mn-6 nm2+12n2m-12n+(2n-2m+4p-6) (ts+ ntw)。
MPI源程序请参见所附光盘。
奇异值分解
A=[aij] 为一个m×n阶矩阵，若存在各列两两正交的m×n的矩阵U、n阶正交方阵V与主对角元素全为非负的对角阵D=diag (d0, d1,…, d(n-1))使A=UDVT，则称d0, d1,…, d(n-1)为A的奇异值。将A写成上式右端的形式，称为对A进行奇异值分解(Singular Value Decomposition)。其中U的各列向量u0, u1,…, u(n-1)为A的左奇异向量，V的各列向量v0 , v1,…, v(n-1)为A的右奇异向量。
矩阵奇异值分解的串行算法
这里我们将介绍对矩阵进行奇异值分解的Henstenes方法。它的基本思想是产生一个正交方阵V，使AV=Q，Q的各列之间也两两正交，若将Q的各列标准化，即使各列向量qi的长度(qi,qi)1/2等于1，这里(s,t)表示向量s, t 的内积。这样，得到：
Q = UD
此处，U满足：UTU =I，D为对角阵：D=diag (d0, d1, …, d(n-1))，非负数di为Q第i列的长度。
由Q=UD可得到：A= QVT =UDVT，即A的奇异值分解。
Henstenes方法应用逐次平面旋转来求Q与V。设A=A0，逐次选择旋转方阵Rk，使Ak的某两列正交化，反复进行旋转变换：
Ak+1 = AkRk (k=0,1,…)
可得矩阵序列{Ak}，适当选择列的正交化顺序，Ak的各列将趋向于两两正交，即Ak趋向于Q，而V=R0,R1,…, Rs，这里s为旋转的总次数。
设Rk使Ak的第p,q两列正交，(p<q)，且Rk =[rij]，取：
rpp = rqq = cosθ rpq = - rqp = sinθ
rii = 1 (i ≠ p, i ≠ q) rij = 0 (i ≠ p, q, j ≠ p,q, i ≠ j )
显然Rk为正交方阵，我们称Rk为Givens旋转矩阵，这样的变换为Givens旋转变换（Givens Rotation）。Ak+1仅在第p,q两列上与Ak不同：
即
为了使Ak+1的 p, q列能够正交，旋转角度θ必须按如下公式选取：假设，若e=0则说明此两列已经正交，取θ=0；否则记 ，由此可求得θ，满足。
被正交化的两列及为一个列对，简记为[p,q]，用上述旋转方法对所有列对[p,q]( p<q)按照i=0,1,…,n-1, j=i+1,…,n-1的顺序正交化一次，称为一轮，一轮计算以后，再做第二轮计算直至A的各列两两正交为止。实际计算中，可给定一个足够小的正数ε，当所有列对的内积e都满足│e│<ε时，计算即可结束。奇异值分解串行算法如下：
算法18.13 矩阵奇异值分解串行算法
输入：矩阵A和矩阵E
输出：矩阵U, 矩阵D和矩阵V T
Begin
(1)while (p>ε)
p=0
for i=1 to n do
for j=i+1 to n do
(1.1 )sum[0]=0 , sum[1]=0 , sum[2]=0
(1.2 )for k=1 to m do
sum[0]= sum[0]+a[k,i]*a[k, j]
sum[1]= sum[1]+ a[k,i]* a[k,i]
sum[2]= sum[2]+ a[k, j]* a[k, j]
end for
(1.3 )if (│sum[0]│>ε) then
(i)aa=2*sum[0]
bb=sum[1]-sum[2]
rr=sqrt(aa*aa+bb*bb)
(ii)if (bb ≥ 0) then
c=sqrt((bb+rr)/(2*rr))
s=aa/(2*rr*c)
else
s=sqrt((rr-bb)/(2*rr))
c=aa/(2*rr*s)
end if
(iii)for k=1 to m do
temp[k]=c*a[k,i]+s*a[k,j]
a[k,j] =(-s)*a[k,i]+c*a[k,j]
end for
(iv)for k=1 to n do
temp1[k]= c*e[k,i]+s*e[k,j]
e[k,j]=(-s)*e[k,i]+c*e[k,j]
end for
(v) for k=1 to m do
a[k,i]= temp[k]
end for
(vi)for k=1 to n do
e[k,i]= temp1[k]
end for
(vii)if (│sum[0]│>p) then p=│sum[0] │ end if
end if
end for
end for
end while
(2)for i=1 to n do
(2.1)sum=0
(2.2)for j=1 to m do
sum=sum+a[j,i]* a[j,i]
end for
(2.3) D[i,i]=sqrt[sum]
end for
(3)for j=1 to n do
for i=1 to m do
U[i,j]=A[i,j]/D[j,j]
end for
end for
(4)VT=MATRIX_TRANSPOSITION(E) /*对E实施算法18.1的矩阵转置*/
End
上述算法的一轮迭代需进行n(n-1)/2次旋转变换，若取一次乘法或一次加法运算时间为一个单位时间，则一次旋转变换要做3次内积运算而消耗6m个单位时间；与此同时，两列元素进行正交计算还需要6m+6n个单位时间，所以奇异值分解的一轮计算时间复杂度为n(n-1)/2*(12m+6n)= n(n-1) (6m+3n)=O(n2m)。
矩阵奇异值分解的并行算法
在并行计算矩阵奇异值分解时，对m×n的矩阵A按行划分为p块(p为处理器数)，每块含有连续的q行向量，这里，第i块包含A的第i×q, …,(i+1)×q-1行向量，其数据元素被分配到第i号处理器上(i=0,1,…,p-1)。E矩阵取阶数为n的单位矩阵I，按同样方式进行数据划分，记，每块含有连续的z行向量。对矩阵A的每一个列向量而言，它被分成p个长度为q的子向量，分布于p个处理器之中，它们协同地对各列向量做正交计算。在对第i列与第j列进行正交计算时，各个处理器首先求其局部存储器中的q维子向量[i,j]、[i,i]、[j,j]的内积，然后通过归约操作的求和运算求得整个m维向量对[i,j]、[i,i]、[j,j]的内积并广播给所有处理器，最后各处理器利用这些内积对其局部存储器中的第i列及第j列q维子向量的元素做正交计算。下述算法18.14按这样的方式对所有列对正交化一次以完成一轮运算，重复进行若干轮运算，直到迭代收敛为止。具体算法框架描述如下：
算法18.14 矩阵奇异值分解并行算法
输入：矩阵A和矩阵E
输出：矩阵U, 矩阵D和矩阵V T
Begin
对所有处理器my_rank(my_rank=0,…, p-1)同时执行如下的算法:
while (not convergence) do
(1)k=0
(2)for i=0 to n-1 do
for j=i+1 to n-1 do
(2.1)sum[0]=0, sum[1]=0, sum[2]=0
(2.1)计算本处理器所存的列子向量a[*,i]、a[*,j]的内积赋给ss[0]
(2.3)计算本处理器所存的列子向量a[*,i]、a[*,j]的内积赋给ss[1]
(2.4)计算本处理器所存的列子向量a[*,i]、a[*,j]的内积赋给ss[2]
(2.5)通过规约操作实现:
(i)计算所有处理器中ss[0]的和赋给 sum[0]
(ii)计算所有处理器中ss[1]的和赋给sum[1]
(iii)计算所有处理器中ss[2]的和赋给sum[2]
(iv)将sum [0]、sum [1]、sum [2] 的值广播到所有处理器中
(2.6)if (│sum(0)│>ε) then /*各个处理器并行地对i和j列的正交化*/
(i)aa=2*sum[0]
(ii)bb=sum[1]-sum[2]
(iii)rr=sqrt(aa*aa+bb*bb)
(iv)if (bb>=0) then
c=sqrt((bb+rr)/(2*rr))
s=aa/(2*rr*c)
else
s=sqrt((rr-bb)/(2*rr))
c=aa/(2*rr*s)
end if
(v)for v=0 to q-1 do
temp[v]=c*a[v, i]+s*a[v, j]
a[v, j] =(-s)*a[v, i]+c*a[v, j]
end for
(vi)for v=0 to z-1 do
temp1[v]=c*e[v, i]+s*e[v, j]
e[v,j] =(-s)*e[v, i]+c*e[v, j]
end for
(vii)for v=0 to q-1 do
a[v,i]=temp[v]
end for
(viii)for v=0 to z-1 do
e[v,i]=temp1[v]
end for
(ix ) k=k+1
end if
end for
end for
end while
End
计算完成后，编号为0的处理器收集各处理器中的计算结果，并进一步计算出矩阵U，D和VT。算法18.14一轮迭代要进行n(n-1)/2次旋转变换,一次旋转变换需要做3次内积运算，若取一次乘法或一次加法运算时间为一个单位时间，则需要6q个单位时间，另外，还要对四列子向量中元素进行正交计算花费6q+6z个单位时间，所以一轮迭代需要的计算时间为n(n-1)(6q+3z)；内积计算需要通信，一轮迭代共做归约操作n(n-1)/2次，每次通信量为3，因而通信时间为。由此得出一轮迭代的并行计算时间为。
MPI源程序请参见所附光盘。
Cholesky分解
对于n阶方阵A=[aij]，若满足A=AT，则称方阵A为对称矩阵(Symmetric Matrix)。对非奇异对称方阵的LU分解有其特殊性，由线性代数知识可知：对于一个对称方阵A，存在一个下三角方阵L，使得A=LLT。由A求得L的过程，称为对A的Cholesky分解。
矩阵Cholesky分解的串行算法
在对矩阵A进行Cholesky分解时,记L=[lij]，当i<j时，lij=0，当i>j时，lij可由下列递推式求得：
aij(1)=aij
aij(k+1)= aij(k)+lik(-ljk)
lij=aij(j)/ljj
当i=j时，lii则可由下列递推式求得：
aii(1)=aii
aii(k+1)= aii(k)-lik2
在串行计算时，可按L的行顺序逐行计算其各个元素，如果取一次乘法和加法运算时间或一次除法运算时间为一个单位时间，则下述算法18.15的计算时间为。
算法18.15 单处理器上矩阵的Cholesky分解的串行算法
输入：矩阵An×n
输出：下三角矩阵Ln×n
Begin
(1)for k=1 to n do
(1.1)a[k, k]=sqrt(a[k, k])
(1.2)for i=k+1 to n do
(i)a[i, k]=a[i, k]/a[k, k]
(ii)for j=k+1 to i do
a[i, j]=a[i, j]-a[i, k]*a[j, k]
end for
end for
end for
(2)for i= 1 to n do
for j= 1 to n do
if (j ≤ i) then l[i,j]= a[i, j] end if
end for
end for
End
矩阵Cholesky分解的并行算法
设A=[aij]，G =[gij]均为n×n阶矩阵，A为实对称正定矩阵，G为上三角矩阵且A=GTG，设则有：
即
令。A′和G′的阶数均减小了1，这样进行n-1次则可以使问题化简到1阶的情况。又注意到A、G都是对称矩阵所以可以考虑节省空间的压缩存储法将A的第k次操作矩阵记作则Cholersky 分解的变换公式推导如下：
算法18.16 矩阵Cholesky分解的并行算法
输入：矩阵An×n,
输出：对应的上对角阵
Begin
对所有处理器my_rank(my_rank=0,…, p-1)同时执行如下的算法:
while (k<n) do
(1)将aij(i=k,j=k,…, n)分配给p个处理器，让每个处理器计算
(2)聚集相应数据，且将第k列对应于第k行更新，将它们广播到各处理器
(3)将aij(i=k+1,…, n, j=i+1,…, n)分配给p个处理器同时计算
(4)聚集(3)中各处理器的数据并广播到各处理器中
(5)k=k+1
(6)
end while
End
很显然，本算法中开平方的次数可减少到n，而乘除次数为，所以时间复杂度为。
MPI源程序请参见所附光盘。
方阵求逆
矩阵求逆(Matrix Inversion)是一常用的矩阵运算。对于一个n×n阶的非奇异方阵A=[aij]，其逆矩阵是指满足A-1A= AA-1=I的n×n阶方阵，其中I为单位方阵。
求方阵的逆的串行算法
这里，我们不妨记A(0)=A，用Gauss-Jordan消去法求A-1的算法是由A(0)出发通过变换得到方阵序列{A(k)}，A(k) ＝ [aij(k)]，(k=0,1,…,n)，每次由A(k-1)到A(k)的变换执行下述的计算：
对于k行的其它诸元素： j=1,2, …,n;jk
除k行、k列外的其它元素： (1i,jn, ik, jk)
对于k列的其它诸元素： i=1,2, …,n;ik
这样计算得到的A(n)就是A-1，矩阵求逆串行算法如算法18.17所示。假定一次乘法和加法运算时间或一次除法运算时间为一个单位时间，则下述矩阵求逆算法18.17的时间复杂度为О(n3)。
算法18.17 单处理器上的矩阵求逆算法
输入：矩阵An×n
输出：矩阵A-1n×n
Begin
for i=1 to n do
(1) a[i,i]=1/a[i,i]
(2)for j=1 to n do
if (j≠i) then a[i,j]=a[i,j]*a[i,i] end if
end for
(3)for k=1 to n do
for j=1 to n do
if ((k≠i and j≠i)) then a[k,j]=a[k,j]-a[k,i]*a[i,j] end if
end for
end for
(4)for k=1 to n do
if (k≠i) then a[k,i]= -a[k,i]*a[i,i] end if
end for
end for
End
方阵求逆的并行算法
矩阵求逆的过程中，依次利用主行i(i=0,1,…,n-1)对其余各行j(j≠i)作初等行变换，由于各行计算之间没有数据相关关系，因此我们对矩阵A按行划分来实现并行计算。考虑到在计算过程中处理器之间的负载均衡，对A采用行交叉划分：设处理器个数为p，矩阵A的阶数为n，，对矩阵A行交叉划分后，编号为i(i=0,1,…,p-1)的处理器存有A的第i, i+p,…, i+(m-1)p行。在计算中，依次将第0,1,…,n-1行作为主行，将其广播给所有处理器，这实际上是各处理器轮流选出主行并广播。发送主行数据的处理器利用主行对其主行之外的m-1行行向量做行变换，其余处理器则利用主行对其m行行向量做行变换。具体算法框架描述如下：
算法18.18 矩阵求逆的并行算法
输入：矩阵An×n,
输出：矩阵A-1n×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为主元素的行号*/
a[i,v]=1/ a[i,v]
(1.2)for k=0 to n-1 do
if (k≠v) then a[i,k]=a[i,k]*a[i,v] end if
end for
(1.3)for k= 0 to n-1 do
f[k]=a[i,k]
end for
(1.4)将变换后的主行元素(存于数组f中)广播到所有处理器中
else /* 主元素不在本处理器中*/
(1.5) v=i*p+j /* v为主元素的行号*/
(1.6) 接收广播来的主行元素存于数组f中
end if
(2)if (my_rank≠j) then /* 主元素不在本处理器中*/
(2.1)for k= 0 to m-1 do /*处理非主行、非主列元素*/
for w= 0 to n-1 do
if (w ≠ v) then a[k,w]=a[k,w]-f[w]*a[k,v] end if
end for
end for
(2.2)for k= 0 to m-1 do /*处理主列元素*/
a[k,v]=-f[v]*a[k,v]
end for
else /*处理主行所在的处理器中的其它元素*/
(2.3)for k= 0 to m-1 do
if ( k ≠ i) then
for w= 0 to n-1 do
if (w ≠ v) then a[k,w]=a[k,w]-f[w]*a[k,v] end if
end for
end if
end for
(2.4)for k= 0 to m-1 do
if ( k ≠ i) then a[k,v]=-f[v]*a[k,v] end if
end for
end if
end for
end for
End
若一次乘法和加法运算或一次除法运算时间为一个单位时间，则算法18.18所需的计算时间为mn2；又由于共有n行数据依次作为主行被广播，其通信时间为n (ts+ ntw)logp，所以该算法并行计算时间为Tp=mn2+n (ts+ ntw)logp。
MPI源程序请参见章末附录。
小结
本章讨论了矩阵转置、矩阵向量乘、矩阵乘法、矩阵的分解及其它一些有关矩阵的数值计算问题，这方面的更详尽的讲解可参见[1]、[2]。关于并行与分布式数值算法，[3]是一本很好的参考书。有关矩阵运算，[4]是一本比较经典的教本。本章所讨论的Cannon乘法原始论文来源于[5]，习题中的Fox乘法来源于[6]，这些算法的改进版本可参见[7]和[8]。
参考文献
陈国良 编著．并行计算——结构·算法·编程．高等教育出版社,1999.10
陈国良 编著．并行算法的设计与分析（修订版）．高等教育出版社，2002.11
Bertsekas D P and Tsitsilklis J N.Parallel and Distributed Computation:Numerical Methods. Prentice-Hall,1989
Golub G H,Loan C V.Matrix Computations.(2ndEd).The Johns Hopkings Univ.Press,1989
Cannon L E.A Cellular Computer to Implement the Kalman Filter Algorithm: Ph.D.thesis. Montana State Univ.,1969
Fox G C,Otto S W,Hey A J G.Matrix Algorithms on Hypercube I:Matrix Multiplication.Parallel Computing,1987,4:17~31
Berntsen J.Communication Effcient Matrix Multiplication on Hypercubes.Parallel Computing,1989,12:335~342
Ho C T,Johnsson S L,Edelman A.Matrix Multiplication on Hypercubes using Full Bandwidth and Constant Storage.Proc.Int’1 91 Conference on Parallel Processing,1997,447~451
附录一. Cannon乘法并行算法的MPI源程序
1. 源程序cannon.cc
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "mpi.h"
#define MAX_PROCESSOR_NUM 12
void MatrixMultiply(int n, double *a, double *b, double *c);
main (int argc, char *argv[])
{
int i,j,k,m,p;
int n, nlocal;
double *a, *b, *c;
int npes, dims[2], periods[2];
int myrank, my2drank, mycoords[2];
int shiftsource, shiftdest, rightrank;
int leftrank, downrank, uprank;
MPI_Status status;
MPI_Comm comm_2d;
MPI_Init(&argc, &argv);
MPI_Comm_size(MPI_COMM_WORLD,
&npes);
MPI_Comm_rank(MPI_COMM_WORLD,
&myrank);
if (argc !=2)
{
if (myrank ==0)
printf("Usage: %s <the dimension of
thematrix>\n", argv[0]);
MPI_Finalize();
exit(0);
}
n=atoi(argv[1]);
dims[0] = sqrt(npes);
dims[1] = npes/dims[0];
if (dims[0] != dims[1])
{
if (myrank ==0 )
printf("The number of processes must be a
perfect square.\n");
MPI_Finalize();
exit(0);
}
periods[0]=periods[1]=1;
MPI_Cart_create(MPI_COMM_WORLD,2,
dims, periods, 0, &comm_2d);
MPI_Comm_rank(comm_2d, &my2drank);
MPI_Cart_coords(comm_2d,my2drank,2,
mycoords);
nlocal = n/dims[0];
a=(double*)malloc(nlocal*nlocal*
sizeof (double));
b= (double *)malloc(nlocal*nlocal*
sizeof (double));
c= (double *)malloc(nlocal*nlocal*
sizeof (double));
srand48((long)myrank);
for ( i=0; i<nlocal*nlocal; i++)
{
a[i] = b[i] =(double)i;
c[i] = 0.0;
}
MPI_Cart_shift(comm_2d, 0, -mycoords[0],
&shiftsource, &shiftdest);
MPI_Sendrecv_replace(a,nlocal*nlocal,
MPI_DOUBLE, shiftdest,1,shiftsource, 1,
comm_2d, &status);
MPI_Cart_shift(comm_2d, 1, -mycoords[1],
&shiftsource, &shiftdest);
MPI_Sendrecv_replace(b, nlocal*nlocal,
MPI_DOUBLE,shiftdest, 1, shiftsource, 1,
comm_2d, &status);
MPI_Cart_shift(comm_2d, 0,-1, &rightrank,
&leftrank);
MPI_Cart_shift(comm_2d, 1, -1, &downrank,
&uprank);
for (i=0; i<dims[0]; i++)
{
MatrixMultiply(nlocal, a, b, c);
MPI_Sendrecv_replace(a, nlocal*nlocal,
MPI_DOUBLE,leftrank,1,rightrank,
1, comm_2d, &status);
MPI_Sendrecv_replace(b, nlocal*nlocal,
MPI_DOUBLE,uprank,1, downrank,
1, comm_2d, &status);
}
MPI_Cart_shift(comm_2d, 0, +mycoords[0],
&shiftsource, &shiftdest);
MPI_Sendrecv_replace(a, nlocal*nlocal,
MPI_DOUBLE, shiftdest, 1, shiftsource,
1, comm_2d, &status);
MPI_Cart_shift(comm_2d, 1, +mycoords[1],
&shiftsource, &shiftdest);
MPI_Sendrecv_replace(b, nlocal*nlocal,
MPI_DOUBLE, shiftdest, 1, shiftsource,
1,comm_2d, &status);
MPI_Comm_free(&comm_2d);
if (myrank ==0)
{
puts("Random Matrix A");
for(i = 0; i < nlocal; i ++)
{
for(j = 0; j < nlocal; j ++)
printf("%9.7f ", a[i*nlocal+j]);
printf("\n");
}
puts("Random Matrix B");
for(i = 0; i < nlocal; i ++)
{
for(j = 0; j < nlocal; j ++)
printf("%9.7f ", b[i*nlocal+j]);
printf("\n");
}
puts("Matrix C = A*B");
for(i = 0; i < nlocal; i ++)
{
for(j = 0; j < nlocal; j ++)
printf("%9.7f ", c[i*nlocal+j]);
printf("\n");
}
}
free(a);
free(b);
free(c);
MPI_Finalize();
}
void MatrixMultiply(int n, double *a, double *b, double *c)
{
int i, j, k;
for (i = 0; i < n; i++)
for (j = 0; j < n; j++)
for (k = 0; k < n; k++)
c[i*n+j]+=a[i*n+k]*b[k*n+j];
} 2. 运行实例
编译：mpicc cannon.cc –o cannon
运行：可以使用命令 mpirun –np ProcessSize cannon ArraySize来运行该Cannon乘法程序，其中ProcessSize是所使用的处理器个数,ArraySize是矩阵的维数,这里分别取为3和4。本实例中使用了
mpirun –np 3 cannon 4
运行结果：
Random Matrix A
0.1708280 0.7499020
0.0963717 0.8704652
Random Matrix B
0.1708280 0.7499020
0.0963717 0.8704652
Matrix C = A*B
0.3999800 1.1517694
1.0546739 1.2320793
说明：该运行实例中，A为2×2的矩阵，B为2×2的矩阵，两者相乘的结果存放于矩阵C中输出。
附录二. 矩阵LU分解并行算法的MPI源程序
1. 源程序 ludep.c
#include "stdio.h"
#include "stdlib.h"
#include "mpi.h"
#define a(x,y) a(x*M+y)
/*A为M*M矩阵*/
#define A(x,y) A(x*M+y)
#define l(x,y) l(x*M+y)
#define u(x,y) u(x*M+y)
#define floatsize sizeof(float)
#define intsize sizeof(int)
int M,N;
int m;
float *A;
int my_rank;
int p;
MPI_Status status;
void fatal(char *message)
{
printf("%s\n",message);
exit(1);
}
void Environment_Finalize(float *a,float *f)
{
free(a);
free(f);
}
int main(int argc, char **argv)
{
int i,j,k,my_rank,group_size;
int i1,i2;
int v,w;
float *a,*f,*l,*u;
FILE *fdA;
MPI_Init(&argc,&argv);
MPI_Comm_size(MPI_COMM_WORLD,
&group_size);
MPI_Comm_rank(MPI_COMM_WORLD,
&my_rank);
p=group_size;
if (my_rank==0)
{
fdA=fopen("dataIn.txt","r");
fscanf(fdA,"%d %d", &M, &N);
if(M != N)
{
puts("The input is error!");
exit(0);
}
A=(float *)malloc(floatsize*M*M);
for(i = 0; i < M; i ++)
for(j = 0; j < M; j ++)
fscanf(fdA, "%f", A+i*M+j);
fclose(fdA);
}
/*0号处理器将M广播给所有处理器*/
MPI_Bcast(&M,1,MPI_INT,0,
MPI_COMM_WORLD);
m=M/p;
if (M%p!=0) m++;
/*分配至各处理器的子矩阵大小为m*M*/
a=(float*)malloc(floatsize*m*M);
/*各处理器为主行元素建立发送和接收缓冲
区*/
f=(float*)malloc(floatsize*M);
/*0号处理器为l和u矩阵分配内存，以分离
出经过变换后的A矩阵中的l和u矩阵
*/
if (my_rank==0)
{
l=(float*)malloc(floatsize*M*M);
u=(float*)malloc(floatsize*M*M);
}
/*0号处理器采用行交叉划分将矩阵A划分为
大小m*M的p块子矩阵，依次发送给1
至p-1号处理器*/
if (a==NULL) fatal("allocate error\n");
if (my_rank==0)
{
for(i=0;i<m;i++)
for(j=0;j<M;j++)
a(i,j)=A((i*p),j);
for(i=0;i<M;i++)
if ((i%p)!=0)
{
i1=i%p;
i2=i/p+1;
MPI_Send(&A(i,0),M,MPI_FLOAT,
i1,i2,MPI_COMM_WORLD);
}
}
else
{
for(i=0;i<m;i++)
MPI_Recv(&a(i,0),M,MPI_FLOAT,
0,i+1,MPI_COMM_WORLD,
&status);
}
for(i=0;i<m;i++)
for(j=0;j<p;j++)
{
/*j号处理器负责广播主行元素*/
if (my_rank==j)
{
v=i*p+j;
for (k=v;k<M;k++)
f(k)=a(i,k);
MPI_Bcast(f,M,MPI_FLOAT,
my_rank,
MPI_COMM_WORLD);
}
else
{
v=i*p+j;
MPI_Bcast(f,M,MPI_FLOAT,
j,
MPI_COMM_WORLD);
}
/*编号小于my_rank的处理器（包
括my_rank本身）利用主行对
其第i+1,…,m-1行数据做行变换*/
if (my_rank<=j)
for(k=i+1;k<m;k++)
{
a(k,v)=a(k,v)/f(v);
for(w=v+1;w<M;w++)
a(k,w)=a(k,w)
-f(w)*a(k,v);
}
/*编号大于my_rank的处理器利用
主行对其第i,…,m-1行数据做行变换*/
if (my_rank>j)
for(k=i;k<m;k++)
{
a(k,v)=a(k,v)/f(v);
for(w=v+1;w<M;w++)
a(k,w)=a(k,w)-f(w)
*a(k,v);
}
}
/*0号处理器从其余各处理器中接收子矩阵a，
得到经过变换的矩阵A*/
if (my_rank==0)
{
for(i=0;i<m;i++)
for(j=0;j<M;j++)
A(i*p,j)=a(i,j);
}
if (my_rank!=0)
{
for(i=0;i<m;i++)
MPI_Send(&a(i,0),M,MPI_FLOAT,
0,i,
MPI_COMM_WORLD);
}
else
{
for(i=1;i<p;i++)
for(j=0;j<m;j++)
{
MPI_Recv(&a(j,0),M,
MPI_FLOAT,i,j,
MPI_COMM_WORLD,
&status);
for(k=0;k<M;k++)
A((j*p+i),k)=a(j,k);
}
}
if (my_rank==0)
{
for(i=0;i<M;i++)
for(j=0;j<M;j++)
u(i,j)=0.0;
for(i=0;i<M;i++)
for(j=0;j<M;j++)
if (i==j)
l(i,j)=1.0;
else
l(i,j)=0.0;
for(i=0;i<M;i++)
for(j=0;j<M;j++)
if (i>j)
l(i,j)=A(i,j);
else
u(i,j)=A(i,j);
printf("Input of file \"dataIn.txt\"\n");
printf("%d\t %d\n",M, N);
for(i=0;i<M;i++)
{
for(j=0;j<N;j++)
printf("%f\t",A(i,j));
printf("\n");
}
printf("\nOutput of LU operation\n");
printf("Matrix L:\n");
for(i=0;i<M;i++)
{
for(j=0;j<M;j++)
printf("%f\t",l(i,j));
printf("\n");
}
printf("Matrix U:\n");
for(i=0;i<M;i++)
{
for(j=0;j<M;j++)
printf("%f\t",u(i,j));
printf("\n");
}
}
MPI_Finalize();
Environment_Finalize(a,f);
return(0);
}
2. 运行实例
编译：mpicc ludep.cc –o ludep
运行：可以使用命令 mpirun –np ProcessSize ludep来运行该矩阵LU分解程序，其中ProcessSize是所使用的处理器个数,这里取为4。本实例中使用了
mpirun –np 4 ludep
运行结果：
Input of file "dataIn.txt"
3 3
2.000000 1.000000 2.000000
0.500000 1.500000 2.000000
2.000000 -0.666667 -0.666667
Output of LU operation
Matrix L:
1.000000 0.000000 0.000000
0.500000 1.000000 0.000000
2.000000 -0.666667 1.000000
Matrix U:
2.000000 1.000000 2.000000
0.000000 1.500000 2.000000
0.000000 0.000000 -0.666667
说明：该运行实例中， A为3×3的矩阵，其元素值存放于文档“dataIn.txt”中，LU分解后得到矩阵L、U作为结果输出。
附录三. 方阵求逆并行算法的MPI源程序
1. 源程序 invert.c
#include "stdio.h"
#include "stdlib.h"
#include "mpi.h"
#define a(x,y) a[x*M+y]
/*A为M*M矩阵*/
#define A(x,y) A[x*M+y]
#define floatsize sizeof(float)
#define intsize sizeof(int)
int M,N;
float *A;
int my_rank;
int p;
MPI_Status status;
void fatal(char *message)
{
printf("%s\n",message);
exit(1);
}
void Environment_Finalize(float *a,float *f)
{
free(a);
free(f);
}
int main(int argc, char **argv)
{
int i,j,k,my_rank,group_size;
int i1,i2;
int v,w;
int m;
float *f;
float *a;
FILE *fdA;
MPI_Init(&argc,&argv);
MPI_Comm_size(MPI_COMM_WORLD,
&group_size);
MPI_Comm_rank(MPI_COMM_WORLD,
&my_rank);
p=group_size;
if (my_rank==0)
{
fdA=fopen("dataIn.txt","r");
fscanf(fdA,"%d %d", &M, &N);
if(M != N)
{
puts("The input is error!");
exit(0);
}
A=(float *)malloc(floatsize*M*M);
for(i = 0; i < M; i ++)
{
for(j = 0; j < M; j ++)
fscanf(fdA, "%f", A+i*M+j);
}
fclose(fdA);
printf("Input of file \"dataIn.txt\"\n");
printf("%d\t%d\n", M, M);
for(i=0;i<M;i++)
{
for(j=0;j<M;j++)
printf("%f\t",A(i,j));
printf("\n");
}
}
/*0号处理器将M广播给所有处理器*/
MPI_Bcast(&M,1,MPI_INT,0,
MPI_COMM_WORLD);
m=M/p;
if (M%p!=0) m++;
/*各处理器为主行元素建立发送和接收缓冲
区*/
f=(float*)malloc(sizeof(float)*M);
/*分配至各处理器的子矩阵大小为m*M*/
a=(float*)malloc(sizeof(float)*m*M);
if (a==NULL||f==NULL)
fatal("allocate error\n");
/*0号处理器采用行交叉划分将矩阵A划分为
m*M的p块子矩阵，依次发送给1至p-1号处理器*/
if (my_rank==0)
{
for(i=0;i<m;i++)
for(j=0;j<M;j++)
a(i,j)=A(i*p,j);
for(i=0;i<M;i++)
if ((i%p)!=0)
{
i1=i%p;
i2=i/p+1;
MPI_Send(&A(i,0),M,
MPI_FLOAT,i1,i2,
MPI_COMM_WORLD);
}
}
else
{
for(i=0;i<m;i++)
MPI_Recv(&a(i,0),M,MPI_FLOAT,
0,i+1,MPI_COMM_WORLD,
&status);
}
for(i=0;i<m;i++)
for(j=0;j<p;j++)
{
/*j号处理器负责广播主行元素*/
if (my_rank==j)
{
v=i*p+j;
a(i,v)=1/a(i,v);
for(k=0;k<M;k++)
if (k!=v)
a(i,k)=a(i,k)*a(i,v);
for(k=0;k<M;k++)
f[k]=a(i,k);
MPI_Bcast(f,M,MPI_FLOAT,
my_rank,
MPI_COMM_WORLD);
}
else
{
v=i*p+j;
MPI_Bcast(f,M,MPI_FLOAT,j
, MPI_COMM_WORLD);
}
/*其余处理器则利用主行对其m行行向量做行变换*/
if (my_rank!=j)
{
for(k=0;k<m;k++)
for(w=0;w<M;w++)
if (w!=v)
a(k,w)=a(k,w)
-f[w]*a(k,v);
for(k=0;k<m;k++)
a(k,v)=-f[v]*
a(k,v);
}
/*发送主行数据的处理器利用主行
对其主行之外的m-1行行向量做行变换*/
if (my_rank==j)
{
for(k=0;k<m;k++)
if (k!=i)
{ for(w=0;w<M;w++)
if (w!=v)
a(k,w)=a(k,w)-
f[w]*a(k,v);
}
for(k=0;k<m;k++)
if (k!=i)
a(k,v)=-f[v]*
a(k,v);
}
}
/*0号处理器从其余各处理器中接收子矩阵a，
得到经过变换的逆矩阵A*/
if (my_rank==0)
{
for(i=0;i<m;i++)
for(j=0;j<M;j++)
A(i*p,j)=a(i,j);
}
if (my_rank!=0)
{
for(i=0;i<m;i++)
for(j=0;j<M;j++)
MPI_Send(&a(i,j),1,
MPI_FLOAT,0,my_rank,
MPI_COMM_WORLD);
}
else
{
for(i=1;i<p;i++)
for(j=0;j<m;j++)
for(k=0;k<M;k++)
{
MPI_Recv(&a(j,k),1,
MPI_FLOAT,i,i,
MPI_COMM_WORLD,
&status);
A((j*p+i),k)=a(j,k);
}
}
if(my_rank==0)
{
printf("\nOutput of Matrix
A'sinversion\n");
for(i=0;i<M;i++)
{
for(j=0;j<M;j++)
printf("%f\t",A(i,j));
printf("\n");
}
}
MPI_Finalize();
Environment_Finalize(a,f);
return(0);
free(A);
}
2. 运行实例
编译：mpicc invert.cc –o invert
运行：可以使用命令 mpirun –np ProcessSize invert来运行该矩阵求逆分解程序，其中ProcessSize是所使用的处理器个数,这里取为3。本实例中使用了
mpirun –np 3 invert
运行结果：
Input of file "dataIn.txt"
3 3
1.000000 -1.000000 1.000000
5.000000 -4.000000 3.000000
2.000000 1.000000 1.000000
Output of Matrix A's inversion
-1.400000 0.400000 0.200000
0.200000 -0.200000 0.400000
2.600000 -0.600000 0.200000
说明：该运行实例中， A为3×3的矩阵，其元素值存放于文档“dataIn.txt”中，求得的逆矩阵A-1作为结果输出。
