1. 6 矩阵运算
矩阵运算是数值计算中最重要的一类运算,特别是在线性代数和数值分析中,它是一种最
基本的运算。本章讨论的矩阵运算包括矩阵转置、矩阵向量相乘、矩阵乘法、矩阵分解以
及方阵求逆等。在讨论并行矩阵算法时分三步进行：①算法描述及其串行算法；②算法的
并行化及其实现算法框架以及简单的算法分析；③算法实现的 MPI 源程序，以利于读者实
践操作。
1.1 矩阵转置
1.1.1 矩阵转置及其串行算法
对于一个 n 阶方阵 A=[a
ij
]，将其每一下三角元素 a
ij
(i>j)沿主对角线与其对称元素 a
ji
互
换就构成了转置矩阵 A
T
。假设一对数据的交换时间为一个单位时间，则下述矩阵转置
(Matrix Transposing)算法 18.1 的运行时间为(n
2
-n)/2=O(n
2
)。
算法 18.1 单处理器上矩阵转置算法
输入：矩阵 A
n×n
输出：矩阵 A
n×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
图 18.1 子块转置
1.1.2 矩阵转置并行算法
此处主要讨论网孔上的块棋盘划分(Block-Checker Board Partitioning，又称为块状划
分)的矩阵转置算法，它对循环棋盘划分(Cyclic-Checker Board Partitioning)也同样适用。另
外，超立方上块棋盘划分的矩阵转置算法可参见文献[1]。
实现矩阵的转置时，若处理器个数为 p，且它们的编号依次是 0,1, …,p-1，则将 n 阶矩
阵 A 分成 p 个大小为 m×m 的子块， 。p 个子块组成一个 的子块阵列。
记其中第 i 行第 j 列的子块为 A
ij
，它含有 A 的第(i-1)m+1 至第 im 行中的第(j-1)m+1 至第 jm
列的所有元素。对每一处理器按行主方式赋以二维下标，记编号为 i 的处理器的二维下标
为(v,u)，其中 ， ，将 A 的子块存入下标为(v,u)表示的对应处理器
中。
这样，转置过程分两步进行：第一步，子块转置，具体过程如图 18 .1 所示；第二步，处
理器内部局部转置。
为了避免对应子块交换数据时处理器发生死锁，可令下三角子块先向与之对应的上三
角子块发送数据，然后从上三角子块接收数据；上三角子块先将数据存放在缓冲区 buffer
中，然后从与之对应的下三角子块接收数据；最后再将缓冲区中的数据发送给下三角子块。
具体并行算法框架描述如下：
算法 18.2 网孔上的矩阵转置算法
输入：矩阵 A
n×n
输出：矩阵 A
n×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
若记 t
s
为发送启动时间, t
w
为单位数据传输时间，t
h
为处理器间的延迟时间，则第一步
由于每个子块有 n
2
/p 个元素，又由于通信过程中为了避免死锁，错开下三角子块与上三角
子块的发送顺序，因此子块的交换时间为 ；第二步，假定一对数据
的交换时间为一个单位时间，则局部转置时间为 。因此所需的并行计算时间
。
MPI 源程序请参见所附光盘。
1.2 矩阵-向量乘法
1.2.1 矩阵-向量乘法及其串行算法
矩阵-向量乘法(Matrix-Vector Multiplication)是将一个 n×n 阶方阵 A=[a
ij
]乘以 n×1 的向量
B=[b
1
,b
2
, …,b
n
]
T
得到一个具有 n 个元素的列向量 C=[c
1
,c
2
, …,c
n
]
T
。假设一次乘法和加法运算
时间为一个单位时间，则下述矩阵向量乘法算法 18.3 的时间复杂度为 O(n
2
)。
算法 18.3 单处理器上矩阵-向量乘算法
输入：A
n×n
，B
n×1
输出：C
n×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 行向量，
，这些行块依次记为 A
0
, A
1
, …, A
p-1
，分别存放在标号为 0,1,…,p-1 的处理器中，
同时将向量 B 广播给所有处理器。各处理器并行地对存于局部数组 a 中的行块 A
i
和向量 B
做乘积操作，具体并行算法框架描述如下：
算法 18.4 行带状划分的矩阵-向量乘并行算法
输入：A
n×n
，B
n×1
输出：C
n×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 并行计算时间 T
p
=n
2
/p，若处理器个数和向量维数相当，则其时间复杂度为 O(n)。
MPI 源程序请参见所附光盘。
1.3 行列划分矩阵乘法
一 个 m×n 阶 矩 阵 A=[a
ij
]乘 以 一个 n×k 的 矩 阵 B=[b
ij
]就 可 以得 到一 个 m×k 的 矩 阵
C=[c
ij
]，它的元素 c
ij
为 A 的第 i 行向量与 B 的第 j 列向量的内积。矩阵相乘的关键是相乘的
两个元素的下标要满足一定的要求(即对准)。为此常采用适当旋转矩阵元素的方法(如后面
将要阐述的 Cannon 乘法)，或采用适当复制矩阵元素的办法(如 DNS 算法)，或采用流水线
的办法使元素下标对准，后两种方法读者可参见文献[1]。
1.3.1 矩阵相乘及其串行算法
由矩阵乘法定义容易给出其串行算法 18.5，若一次乘法和加法运算时间为一个单位时
间，则显然其时间复杂度为 O(mnk) 。
算法 18.5 单处理器上矩阵相乘算法
输入：A
m×n
，B
n×k
输出：C
m×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 行向量，这些行块依次记为 A
0
,A
1
, …,A
p-1
，分别存放在标号为 0,1,…, p-1 的处理器中。
对矩阵 B 按列划分为 P 块，记 ，每块含有连续的 v 列向量，这些列块依次记为
B
0
,B
1
, …,B
p-1
，分别存放在标号 0,1, …, p-1 为的处理器中。将结果矩阵 C 也相应地同时进行
行、列划分，得到 p×p 个大小为 u×v 的子矩阵，记第 i 行第 j 列的子矩阵为 C
ij
，显然有 C
ij
=
A
i
×B
j
，其中，A
i
大小为 u×n，B
j
大小为 n×v。
图 18.2 矩阵相乘并行算法中的数据交换
开始，各处理器的存储内容如图 18 .2 (a)所示。此时各处理器并行计算 C
ii
= A
i
×B
j
其
中 i=0,1,…,p-1，此后第 i 号处理器将其所存储的 B 的列块送至第 i-1 号处理器（第 0 号处理
器将 B 的列块送至第 p-1 号处理器中，形成循环传送），各处理器中的存储内容如 图
18 .2 (b)所示。它们再次并行计算 C
ij
= A
i
×B
j
，这里 j=(i+1)modp。B 的列块在各处理器中以
这样的方式循环传送 p-1 次并做 p 次子矩阵相乘运算，就生成了矩阵 C 的所有子矩阵。编
号为 i 的处理器的内部存储器存有子矩阵 C
i0
,C
i1
,…,C
i(p-1)
。为了避免在通信过程中发生死锁，
奇数号及偶数号处理器的收发顺序被错开，使偶数号处理器先发送后接收；而奇数号处理
器先将 B 的列块存于缓冲区 buffer 中，然后接收编号在其后面的处理器所发送的 B 的列块，
最后再将缓冲区中原矩阵 B 的列块发送给编号在其前面的处理器，具体并行算法框架描述
如下：
算法 18.6 矩阵并行分块乘法算法
输入：A
m×n
, B
n×k
,
输出：C
m×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 中做备份
