并 行 计 算中国科学技术大学计算机科学与技术系国家高性能计算中心 (合肥 )
2003年 9月第三篇 并行数值算法第八章 基本通讯操作第九章 稠密矩阵运算第十章 线性方程组的求解第十一章 快速傅里叶变换第九章 稠密矩阵运算
9.1 矩阵的划分
9.2 矩阵转臵
9.3 矩阵 -向量乘法
9.4 矩阵乘法
9.1 矩阵的划分
9.1.1 带状划分
9.1.2 棋盘划分国家高性能计算中心（合肥） 52009-7-24
带状划分
16× 16阶矩阵,p=4
列块带状划分 行循环带状划分
P P P P
0
4
8
1 2
1
5
9
1 3
2
6
1 0
1 4
3
7
1 1
1 5
P
P
P
P
( b )
( a )
图 9,1
0 1 2 3
3
2
1
0
0 1 2 3 4 5 6 7 8 9 1 0 1 1 1 2 1 3 1 4 1 5
国家高性能计算中心（合肥） 62009-7-24
带状划分
示例,p＝ 3,27× 27矩阵的 3种带状划分
9.1 矩阵的划分
9.1.1 带状划分
9.1.2 棋盘划分国家高性能计算中心（合肥） 82009-7-24
棋盘划分
8× 8阶矩阵,p=16
块棋盘划分 循环棋盘划分
( b )( a )
图 9,2
6
1 0
1 4
7
1 1
1 5
4
8
1 2
5
9
1 3
6
1 0
1 4
( 0,0 )
( 1,0 )
( 2,0 )
( 3,0 )
( 5,0 )
( 4,0 )
( 7,0 )
( 6,0 )
( 0,1 )
( 1,1 )
( 2,1 )
( 3,1 )
( 5,1 )
( 4,1 )
( 7,1 )
( 6,1 )
( 0,2 )
( 1,2 )
( 2,2 )
( 3,2 )
( 5,2 )
( 4,2 )
( 7,2 )
( 6,2 )
( 0,3 )
( 1,3 )
( 0,4 )
( 1,4 )
( 2,3 )
( 3,3 )
( 2,4 )
( 3,4 )
( 5,3 )
( 4,3 )
( 5,4 )
( 4,4 )
( 7,3 )
( 6,3 )
( 7,4 )
( 6,4 )
( 0,5 )
( 1,5 )
( 0,6 )
( 1,6 )
( 2,5 )
( 3,5 )
( 2,6 )
( 3,6 )
( 5,5 )
( 4,5 )
( 5,6 )
( 4,6 )
( 7,5 )
( 6,5 )
( 7,6 )
( 6,6 )
( 0,7 )
( 1,7 )
( 2,7 )
( 3,7 )
( 5,7 )
( 4,7 )
( 7,7 )
( 6,7 )
( 0,0 )
( 2,0 )
( 3,0 )
( 5,0 )
( 4,0 )
( 7,0 )
( 6,0 )
( 1,0 )
( 6,1 )
( 1,1 )
( 0,4 )
( 1,4 )
( 2,1 )
( 3,1 )
( 2,4 )
( 3,4 )
( 5,1 )
( 4,1 )
( 5,4 )
( 4,4 )
( 7,1 )( 7,4 )
( 6,4 )
( 0,1 )
( 1,2 )
( 0,5 )
( 1,5 )
( 2,2 )
( 3,2 )
( 2,5 )
( 3,5 )
( 5,2 )
( 4,2 )
( 5,5 )
( 4,5 )
( 7,2 )
( 6,2 )
( 7,5 )
( 6,5 )
( 0,2 )
( 1,3 )
( 0,6 )
( 1,6 )
( 2,3 )( 2,6 )
( 3,6 )
( 5,3 )
( 4,3 )
( 5,6 )
( 4,6 )
( 7,3 )
( 6,3 )
( 7,6 )
( 6,6 )
( 0,3 )
( 3,3 )
( 0,7 )
( 1,7 )
( 2,7 )
( 3,7 )
( 5,7 )
( 4,7 )
( 7,7 )
( 6,7 )
3
7
1 1
1 5
0
4
8
1 2
1
5
9
1 3
P P P P P P P P
2 3 0 1 2
P P P P P P P P
P P P P P P P P
P P P P P P P P
国家高性能计算中心（合肥） 92009-7-24
棋盘划分
示例,p＝ 4,16× 16矩阵的 3种棋盘划分第九章 稠密矩阵运算
9.1 矩阵的划分
9.2 矩阵转臵
9.3 矩阵 -向量乘法
9.4 矩阵乘法
9.2 矩阵转臵
9.2.1 棋盘划分的矩阵转臵
9.2.2 带状划分的矩阵转臵国家高性能计算中心（合肥） 122009-7-24
棋盘划分的矩阵转置
网孔连接
情形 1,p=n2。
通讯步 转臵后
P
( 3,0 )
( 1,0 ) ( 1,2 ) ( 1,3 )
( 2,0 ) ( 2,1 ) ( 2,3 )
P
( 3,1 )
P
( 3,2 )
P P P P
( b )( a )
图 9,3
( 0,0 )
( 0,1 )
( 0,2 )
( 0,3 )
( 1,0 )
( 1,1 )
( 1,2 )
( 1,3 )
( 2,0 )
( 2,1 )
( 2,2 )
( 2,3 )
( 3,0 )
( 3,1 )
( 3,2 )
( 3,3 )
( 0,1 ) ( 0,2 ) ( 0,3 )
P
( 0,0 )
( 1,1 )
( 2,2 )
( 3,3 )
3
7
1 1
1 5
2
6
1 0
1 4
1
5
9
1 31 2
0
4
8
3
7
1 1
1 5
2
6
1 4
1 0
1
9
1 3
5
P P P P P P PP
P P P P P P PP
P P P P P P PP
P
1 2
4
8
0
国家高性能计算中心（合肥） 132009-7-24
棋盘划分的矩阵转置
情形 2,p<n2。
- 划分,
- 算法,① 按 mesh连接进行块转臵 (不同处理器间 )
② 进行块内转臵 (同一处理器内 )
通讯步 转臵后
P P P
P P P
P P P
P P P
( b )
( a )
图 9,4
P P P P
P P P P
P P P P
P P P P
0
4
8
1 2
1
5
9
1 3
2
6
1 0
1 4
3
7
1 1
1 5
( 0,0 )
( 1,0 )
( 0,2 )
( 1,2 )
( 1,4 )
( 0,4 )
( 1,6 )
( 0,6 )
( 0,1 )
( 1,1 )
( 0,3 )
( 1,3 )
( 1,5 )
( 0,5 )
( 1,7 )
( 0,7 )
( 2,0 )
( 3,0 )
( 2,2 )
( 3,2 )
( 3,4 )
( 2,4 )
( 3,6 )
( 2,6 )
( 2,1 )
( 3,1 )
( 2,3 )
( 3,3 )
( 3,5 )
( 2,5 )
( 3,7 )
( 2,7 )
( 4,0 )
( 5,0 )
( 4,2 )
( 5,2 )
( 5,4 )
( 4,4 )
( 5,6 )
( 4,6 )
( 4,1 )
( 5,1 )
( 4,3 )
( 5,3 )
( 5,5 )
( 4,5 )
( 5,7 )
( 4,7 )
( 6,0 )
( 7,0 )
( 6,2 )
( 7,2 )
( 7,4 )
( 6,4 )
( 7,6 )
( 6,6 )
( 6,1 )
( 7,1 )
( 6,3 )
( 7,3 )
( 7,5 )
( 6,5 )
( 7,7 )
( 6,7 )
4
8
1 2
0
P
P
1
5
9
1 3
P
2
6
1 0
1 4
P
3
7
1 1
1 5
( 1,0 )
( 2,0 )
( 3,0 )
( 5,0 )
( 4,0 )
( 7,0 )
( 6,0 )
( 0,0 ) ( 0,1 )
( 1,1 )
( 2,1 )
( 3,1 )
( 5,1 )
( 4,1 )
( 7,1 )
( 6,1 )
( 0,2 )
( 1,2 )
( 2,2 )
( 3,2 )
( 5,2 )
( 4,2 )
( 7,2 )
( 6,2 )
( 0,3 )
( 1,3 )
( 2,3 )
( 3,3 )
( 5,3 )
( 4,3 )
( 7,3 )
( 6,3 )
( 0,4 )
( 1,4 )
( 2,4 )
( 3,4 )
( 5,4 )
( 4,4 )
( 7,4 )
( 6,4 )
( 0,5 )
( 1,5 )
( 2,5 )
( 3,5 )
( 5,5 )
( 4,5 )
( 7,5 )
( 6,5 )
( 0,6 )
( 1,6 )
( 2,6 )
( 3,6 )
( 5,6 )
( 4,6 )
( 7,6 )
( 6,6 )
( 0,7 )
( 1,7 )
( 2,7 )
( 3,7 )
( 5,7 )
( 4,7 )
( 7,7 )
( 6,7 )
子 块划分成p 个大小为A pnpnnn
通讯?)/(2// 2 pnttp ws?
计算pn2// 2
运行时间pntptpnT wsp /222 22
国家高性能计算中心（合肥） 142009-7-24
棋盘划分的矩阵转置
超立方连接
划分,
算法,
①
② 对 Aij递归应用 ① 进行转臵，直至分块矩阵的元素处于同一处理器；
③ 进行同一处理器的内部转臵。
运行时间,
子 块划分成p 个大小为A pnpnnn
p
p
ntt
p
n
p
p
ntt
p
np
p
ntt
p
nT
ws
wswsp
l o g)
2
l o g)2
2
//l o g)2
2
22
2222
（
递归步：，（选路：，内部转置（
2212 21112221 1211 AA AAAA AAA 转置为＝将国家高性能计算中心（合肥） 152009-7-24
棋盘划分的矩阵转置
超立方连接,示例
( 0,5 )
( 1,5 )
P P P
( 2,5 )
( 3,5 )
P
P P P
( 5,5 )
( 4,5 )
P
P P P
( 7,5 )
( 6,5 )
P
( 0,0 ) ( 0,1 )
( 1,0 ) ( 1,1 )
( 0,2 ) ( 0,3 )
( 1,2 ) ( 1,3 )
P
( 0,4 ) ( 0,5 )
( 1,4 ) ( 1,5 )
P
( 0,6 ) ( 0,7 )
( 1,6 ) ( 1,7 )
P
( 2,0 ) ( 2,1 )
( 3,0 ) ( 3,1 )
P
( 2,2 ) ( 2,3 )
( 3,2 ) ( 3,3 )
P
( 2,4 ) ( 2,5 )
( 3,4 ) ( 3,5 )
P
( 2,6 ) ( 2,7 )
( 3,6 ) ( 3,7 )
( 5,0 ) ( 5,1 )
( 4,0 ) ( 4,1 )
( 5,2 ) ( 5,3 )
( 4,2 ) ( 4,3 )
P
( 5,4 ) ( 5,5 )
( 4,4 ) ( 4,5 )
P
( 5,6 ) ( 5,7 )
( 4,6 ) ( 4,7 )
P
( 7,0 ) ( 7,1 )
( 6,0 ) ( 6,1 )
P
( 7,2 ) ( 7,3 )
( 6,2 ) ( 6,3 )
P
( 7,4 ) ( 7,5 )
( 6,4 ) ( 6,5 )
P
( 7,6 ) ( 7,7 )
( 6,6 ) ( 6,7 )
( b )( a )
图 9,6
0
4
8
1 2
1
5
9
2
6
1 0
1 4
3
7
1 1
1 5
0
8
4
1 2
1
9
5
1 3
6
2
1 0
1 4
7
3
1 1
1 5
( 0,0 )
( 1,0 )
( 2,0 )
( 3,0 )
( 5,0 )
( 4,0 )
( 7,0 )
( 6,0 )
( 0,1 )
( 1,1 )
( 2,1 )
( 3,1 )
( 5,1 )
( 4,1 )
( 7,1 )
( 6,1 )
P P P P P P P P
( 0,2 )
( 1,2 )
( 2,2 )
( 3,2 )
( 5,2 )
( 4,2 )
( 7,2 )
( 6,2 )
( 0,3 )
( 1,3 )
( 2,3 )
( 3,3 )
( 5,3 )
( 4,3 )
1 3
( 7,3 )
( 6,3 )
( 0,4 )
( 1,4 )
( 2,4 )
( 3,4 )
( 5,4 )
( 4,4 )
( 7,4 )
( 6,4 )
( 0,6 )
( 1,6 )
( 2,6 )
( 3,6 )
( 5,6 )
( 4,6 )
( 7,6 )
( 6,6 )
( 0,7 )
( 1,7 )
( 2,7 )
( 3,7 )
( 5,7 )
( 4,7 )
( 7,7 )
( 6,7 )
9.2 矩阵转臵
9.2.1 棋盘划分的矩阵转臵
9.2.2 带状划分的矩阵转臵国家高性能计算中心（合肥） 172009-7-24
带状划分的矩阵转置
划分,An× n分成 p个 (n/p)× n大小的带
算法,
① Pi有 p-1个 (n/p)× (n/p)大小子块发送到另外 p-1个处理器中 ;
② 每个处理器本地交换相应的元素
P
P
P
P
n
图 9,7
0
1
2
3
第九章 稠密矩阵运算
9.1 矩阵的划分
9.2 矩阵转臵
9.3 矩阵 -向量乘法
9.4 矩阵乘法
9.3 矩阵 -向量乘法
9.3.1 带状划分的矩阵 -向量乘法
9.3.2 棋盘划分的矩阵 -向量乘法国家高性能计算中心（合肥） 202009-7-24
带状划分的矩阵 -向量乘法
划分 (行带状划分 ),Pi存放 xi和 ai,0,ai,1,…,ai,n-1,并输出 yi
算法,对 p=n情形
① 每个 Pi向其他处理器播送 xi(多到多播送 )；
② 每个 Pi计算；
注,对 p<n情形,算法中 Pi要播送 X中相应的 n/p个分量
(1)超立方连接的计算时间
(2)网孔连接的计算时间充分大时项是多到多的播送时间后项是乘法时间，前
pntpt
p
n
pt
p
npt
p
nT
ws
wsp
//l o g
21//)1(l o g
2
2
充分大时项是多到多的播送时间后项是乘法时间，前
pntpt
p
n
pt
p
ntp
p
nT
ws
wsp
//)1(2
21//)1()1(2
2
2
国家高性能计算中心（合肥） 212009-7-24
带状划分的矩阵 -向量乘法
示例矩 阵 A
0
1
p - 1
n / p
向 量 x 处 理 器
0
1
p - 1
n
( b )( a )
0
0
0
0
0
1
1
1
1
1 p - 1
矩 阵 A
0
1
p - 1
向 量 y
( d )( c )
图 9,8
P
0
P
1
P
p - 1
P
0
P
1
P
p - 1
P
0
P
1
P
p - 1
P
0
P
1
P
p - 1
p - 1
p - 1
p - 1
p - 1
9.3 矩阵 -向量乘法
9.3.1 带状划分的矩阵 -向量乘法
9.3.2 棋盘划分的矩阵 -向量乘法国家高性能计算中心（合肥） 232009-7-24
棋盘划分的矩阵 -向量乘法
划分 (块棋盘划分 ),Pij存放 ai,j,xi臵入 Pi,i中
算法,对 p=n2情形
① 每个 Pi,i向 Pj,i播送 xi(一到多播送 )；
② 按行方向进行乘 -加与积累运算，最后一列 Pi,n-1收集的结果为 yi；
注,对 p<n2情形,p个处理器排成 的二维网孔，
算法中 Pi,i向 Pj,i播送 X中相应的 个分量
(1)网孔连接的计算时间 Tp(CT),
.X中相应分量臵入 Pi,i的通讯时间,
.按列一到多播送时间,
.按行单点积累的时间,ptpt
p
npt
p
nT
hwsp 3loglog
2
pp?
pn/
pttpnt hws
)1(lo g)( ptptpnt hws
)1(lo g)( ptptpnt hws
国家高性能计算中心（合肥） 242009-7-24
棋盘划分的矩阵 -向量乘法
示例矩 阵 A 向 量 x
n
( d )( c )
图 9,9
P
0
P
1
P
n /
矩 阵 A
向 量
y
P
1
P
0
√ p - 1
P
( b )( a )
√ pn /
√ p - 1
√ p
P
p - 1
P
p - 1
P
P
√ p
2 √ p
P
P
√ p
2 √ p
国家高性能计算中心（合肥） 252009-7-24
带状与棋盘划分比较
以网孔链接为例
网孔上带状划分的运行时间
网孔上棋盘划分的运行时间
棋盘划分要比带状划分快。
)5.9()1(2
2
wsp ntptp
nT
)6.9(3loglog
2
ptptpnptpnT hwsp
第九章 稠密矩阵运算
9.1 矩阵的划分
9.2 矩阵转臵
9.3 矩阵 -向量乘法
9.4 矩阵乘法
9.4 矩阵乘法
9.4.1 简单并行分块乘法
9.4.2 Cannon乘法
9.4.3 Fox乘法
9.4.4 Systolic乘法
9.4.5 DNS乘法国家高性能计算中心（合肥） 282009-7-24
矩阵乘法符号及定义
1,11,10,1
1,11,10,1
1,01,00,0
1,11,10,1
1,11,10,1
1,01,00,0
1,11,10,1
1,11,10,1
1,01,00,0
,)()()(
nnnn
n
n
nnnn
n
n
nnnn
n
n
nnijnnijnnij
bbb
bbb
bbb
aaa
aaa
aaa
ccc
ccc
ccc
BACcCbBaA
设
j
i
A BC
1
0
n
k
kjikij bac
A中元素的第 1下标与 B中元素的第 2下标相一致（对准）
国家高性能计算中心（合肥） 292009-7-24
矩阵乘法并行实现方法
计算结构：二维阵列
空间对准 (元素已加载到阵列中）
Cannon’s,Fox’s,DNS
时间对准 (元素未加载到阵列中）
Systolic A0,0
B0,0
A1,0
B1,0
A2,0
B2,0
A3,0
B3,0
A0,1
B0,1
A1,1
B1,1
A2,1
B2,1
A3,1
B3,1
A0,2
B0,2
A1,2
B1,2
A2,2
B2,2
A3,2
B3,2
A0,3
B0,3
A1,3
B1,3
A2,3
B2,3
A3,3
B3,3
国家高性能计算中心（合肥） 302009-7-24
简单并行分块乘法
分块,A,B和 C分成 的方块阵 Ai,j,Bi,j和 Ci,j,大小均为
p个处理器编号为,Pi,j存放 Ai,j,Bi,j和 Ci,j。
算法,
① 通讯,每行处理器进行 A矩阵块的多到多播送 (得到 Ai,k,k=0~ )
每列处理器进行 B矩阵块的多到多播送 (得到 Bk,j,k=0~ )
② 乘 -加运算,Pi,j做
运行时间
(1)超立方连接：
① 的时间
② 的时间
ppp pnpn?
),.,,,,.,,,( 1,11,00,0 ppp PPP
1－p
10pk kjikij BAC
))1(lo g(2 21 ppntptt ws
pnpt pn /)( 332
1－p
国家高性能计算中心（合肥） 312009-7-24
简单并行分块乘法
运行时间
(1)超立方连接：
(2)二维环绕网孔连接：
① 的时间：
② 的时间 t2=n3/p
注
(1)本算法的缺点是对处理器的存储要求过大每个处理器有 个块，每块大小为 n2/p,
所以需要,p个处理器共需要,
是串行算法的 倍
(2)p=n2时,t(n)=O(n),c(n)=O(n3)
p
ntptpt
p
ntt
wsws
22
1 22)1((2 ）
p
ntpt
p
nT
wsp
23 22
p2
)/( 2 pnO )( 2 pnO
p
9.4 矩阵乘法
9.4.1 简单并行分块乘法
9.4.2 Cannon乘法
9.4.3 Fox乘法
9.4.4 Systolic乘法
9.4.5 DNS乘法国家高性能计算中心（合肥） 332009-7-24
Cannon乘法
分块,A,B和 C分成 的方块阵 Ai,j,Bi,j和 Ci,j,大小均为 p个处理器编号为,Pi,j存放
Ai,j,Bi,j和 Ci,j(n > > p)
ppp
pnpn? ),.,,,,.,,,( 1,11,00,0 ppp PPP
P0,0
P1,0
P2,0
P3,0
P0,1
P1,1
P2,1
P3,1
P0,2
P1,2
P2,2
P3,2
P0,3
P1,3
P2,3
P3,3
p
n
n
p
国家高性能计算中心（合肥） 342009-7-24
Cannon乘法
算法原理 (非形式描述 )
① 所有块 Ai,j(0≤i,j≤ ) 向左循环移动 i步 (按行移位 )；
所有块 Bi,j(0≤i,j≤ ) 向上循环移动 j步 (按列移位 )；
② 所有处理器 Pi,j做执行 Ai,j和 Bi,j的乘 -加运算 ;
③ A的每个块向左循环移动一步 ;
B的每个块向上循环移动一步 ;
④ 转 ② 执行 次 ;
1－p
1－p
1－p
国家高性能计算中心（合肥） 352009-7-24
Cannon乘法
示例,A4× 4,B4× 4,p=16
A0,0
A1,0
A2,0
A3,0
A0,1
A1,1
A2,1
A3,1
A0,2
A1,2
A2,2
A3,2
A0,3
A1,3
A2,3
A3,3
B0,0
B1,0
B2,0
B3,0
B0,1
B1,1
B2,1
B3,1
B0,2
B1,2
B2,2
B3,2
B0,3
B1,3
B2,3
B3,3
Initial alignment of A Initial alignment of B
国家高性能计算中心（合肥） 362009-7-24
Cannon乘法
示例,A4× 4,B4× 4,p=16
A and B after initial alignment and shifts after every step
A0,0
B0,0
A1,1
B1,0
A2,2
B2,0
A3,3
B3,0
A0,1
B1,1
A1,2
B2,1
A2,3
B3,1
A3,0
B0,1
A0,2
B2,2
A1,3
B3,2
A2,0
B0,2
A3,1
B1,2
A0,3
B3,3
A1,0
B0,3
A2,1
B1,3
A3,2
B2,3
国家高性能计算中心（合肥） 372009-7-24
Cannon乘法
示例,A4× 4,B4× 4,p=16
After first shift
A0,1
B1,0
A1,2
B2,0
A2,3
B3,0
A3,0
B0,0
A0,2
B2,1
A1,3
B3,1
A2,0
B0,1
A3,1
B3,1
A0,3
B3,2
A1,0
B0,2
A2,1
B1,2
A3,2
B2,2
A0,0
B0,3
A1,1
B1,3
A2,2
B2,3
A3,3
B3,3
After second shift
A0,2
B2,0
A1,3
B3,0
A2,0
B0,0
A3,1
B1,0
A0,3
B3,1
A1,0
B0,1
A2,1
B1,1
A3,2
B2,1
A0,0
B0,2
A1,1
B1,2
A2,2
B2,2
A3,3
B3,2
A0,1
B1,3
A1,2
B2,3
A2,3
B3,3
A3,0
B0,3
After third shift
A0,3
B3,0
A1,0
B0,0
A2,1
B1,0
A3,2
B2,0
A0,0
B0,1
A1,1
B1,1
A2,2
B2,1
A3,3
B3,1
A0,1
B1,2
A1,2
B2,2
A2,3
B3,2
A3,0
B0,2
A0,2
B2,3
A1,3
B3,3
A2,0
B0,3
A3,1
B1,3
国家高性能计算中心（合肥） 382009-7-24
Cannon乘法
算法描述,Cannon分块乘法算法
//输入,An× n,Bn× n; 输出,Cn× n
Begin
(1)for k=0 to do
for all Pi,j par-do
(i) if i>k then
Ai,j? Ai,(j+1)mod
endif
(ii)if j>k then
Bi,j? B(i+1)mod,j
endif
endfor
endfor
(2)for all Pi,j par-do Ci,j=0 endfor
(3)for k=0 to do
for all Pi,j par-do
(i) Ci,j=Ci,j+Ai,jBi,j
(ii) Ai,j? Ai,(j+1)mod
(iii)Bi,j? B(i+1)mod,j
endfor
endfor
End
时间分析：
)/(
))/(()1()(
)(
3
3
321
pnO
pnpOOpO
TTTnT p
1－p
p
p
p
p
1－p
9.4 矩阵乘法
9.4.1 简单并行分块乘法
9.4.2 Cannon乘法
9.4.3 Fox乘法
9.4.4 Systolic乘法
9.4.5 DNS乘法国家高性能计算中心（合肥） 402009-7-24
Fox乘法
分块,同 Cannon分块算法
算法原理
① Ai,i向所在行的其他处理器进行一到多播送；
② 各处理器将收到的 A块与原有的 B块进行乘 -加运算 ;
③ B块向上循环移动一步 ;
④ 如果 Ai,j是上次第 i行播送的块，本次选择 向所在行的其他处理器进行一到多播送 ;
⑤ 转 ② 执行 次 ;1?p
pjiA m od)1(,?
A0,0
B0,0
A1,0
B1,0
A2,0
B2,0
A3,0
B3,0
A0,1
B0,1
A1,1
B1,1
A2,1
B2,1
A3,1
B3,1
A0,2
B0,2
A1,2
B1,2
A2,2
B2,2
A3,2
B3,2
A0,3
B0,3
A1,3
B1,3
A2,3
B2,3
A3,3
B3,3
国家高性能计算中心（合肥） 412009-7-24
Fox乘法
示例,A4× 4,B4× 4,p=16
(a) (b)
A0,0
B0,0
B1,0
B2,0
B3,0
B0,1
A1,1
B1,1
B2,1
B3,1
B0,2
B1,2
A2,2
B2,2
B3,2
B0,3
B1,3
B2,3
A3,3
B3,3
B1,0
B2,0
B3,0
A0,1
B1,1
B2,1
B3,1
B0,1
B1,2
B3,2
B0,2
B1,3
B2,3
B0,3
A1,2
B2,2
A2,3
B3,3
A3,0
B0,0
国家高性能计算中心（合肥） 422009-7-24
Fox乘法
示例,A4× 4,B4× 4,p=16
(c) (d)
B2,0
B3,0
B2,1
B3,1
B0,1
B3,2
B0,2
B1,2
B2,3
B0,3
B1,3
B3,0
B1,0
B3,1
B0,1
B2,1
B3,2
B1,2
B0,3
B2,3
B0,2
B1,3
B2,0
A0,2
B2,2
A1,3
B3,3
A2,0
B0,0
B1,0
A3,1
B1,1
A0,3
B3,3
A1,0
B0,2
A2,1
B1,1
A3,2
B2,2
9.4 矩阵乘法
9.4.1 简单并行分块乘法
9.4.2 Cannon乘法
9.4.3 Fox乘法
9.4.4 Systolic乘法
9.4.5 DNS乘法国家高性能计算中心（合肥） 442009-7-24
Systolic乘法
a1,4
b4,1
b3,1
b2,1
b2,2
b4,2
b3,2
b2,3
b3,3
b4,3
b2,4
b3,4
b4,4
a1,3a1,1 a1,2
a2,4a2,1 a2,2 a2,3
a3,1 a3,2 a3,3 a3,4
b1,1
b1,2
b1,3
b1,4
Step 1
P1,1
c1,1
P1,2
c1,2
P1,3
c1,3
P1,4
c1,4
P2,1
c2,1
P2,2
c2,2
P2,3
c2,3
P2,4
c2,4
P3,1
c3,1
P3,2
c3,2
P3,3
c3,3
P3,4
c3,4
国家高性能计算中心（合肥） 452009-7-24
Systolic乘法
c1,1 c1,2 c1,3 c1,4
c2,1 c2,2 c2,3 c2,4
c3,1 c3,2 c3,3 c3,4
b3,1
b2,1
b2,2
b4,2
b3,2
b2,3
b3,3
b4,3
b2,4
b3,4
b4,4
a1,3a1,1 a1,2
a2,4a2,1 a2,2 a2,3
a3,1 a3,2 a3,3 a3,4
b1,1
b1,2
b1,3
b1,4
a1,4b4,1
+
Step 2
国家高性能计算中心（合肥） 462009-7-24
Systolic乘法
c1,1 c1,2 c1,3 c1,4
c2,1 c2,2 c2,3 c2,4
c3,1 c3,2 c3,3 c3,4
b2,1
b2,2
b3,2
b2,3
b3,3
b4,3
b2,4
b3,4
b4,4
a1,1 a1,2
a2,1 a2,2 a2,3
a3,1 a3,2 a3,3 a3,4
b1,1
b1,2
b1,3
b1,4
a1,3b3,1
+
a1,4b4,2
+
a2,4b4,1
+
Step 3
国家高性能计算中心（合肥） 472009-7-24
Systolic乘法
c1,1 c1,2 c1,3 c1,4
c2,1 c2,2 c2,3 c2,4
c3,1 c3,2 c3,3 c3,4
b2,2
b2,3
b3,3
b2,4
b3,4
b4,4
a1,1
a2,1 a2,2
a3,1 a3,2 a3,3
b1,1
b1,2
b1,3
b1,4
a1,2b2,1
+
a1,3b3,2
+
a2,3b3,1
+
a1,4b4,3
+
a3,4b4,1
+
a2,4b4,2
+
Step 4
国家高性能计算中心（合肥） 482009-7-24
Systolic乘法
c1,1 c1,2 c1,3 c1,4
c2,1 c2,2 c2,3 c2,4
c3,1 c3,2 c3,3 c3,4
b2,3
b2,4
b3,4
a2,1
a3,1 a3,2
b1,2
b1,3
b1,4
a1,1b1,1
+
a1,2b2,2
+
a2,2b2,1
+
a1,3b3,3
+
a3,3b3,1
+
a2,3b3,2
+
a1,4b4,4
+
a2,4b4,3
+
a3,4b4,2
+
Step 5
国家高性能计算中心（合肥） 492009-7-24
Systolic乘法
c1,1 c1,2 c1,3 c1,4
c2,1 c2,2 c2,3 c2,4
c3,1 c3,2 c3,3 c3,4
b2,4
a3,1
b1,3
b1,4
a1,1b1,2
+
a2,1b1,1
+
a1,2b2,3
+
a3,2b2,1
+
a2,2b2,2
+
a1,3b3,4
+
a2,3b3,3
+
a3,3b3,2
+
a2,4b4,4
+
a3,4b4,3
+
Step 6
国家高性能计算中心（合肥） 502009-7-24
Systolic乘法
c1,1 c1,2 c1,3 c1,4
c2,1 c2,2 c2,3 c2,4
c3,1 c3,2 c3,3 c3,4
b1,4
a1,1b1,3
+
a3,1b1,1
+
a2,1b1,2
+
a1,2b2,4
+
a2,2b2,3
+
a3,2b3,2
+
a2,3b3,4
+
a3,3b3,3
+
a3,4b4,4
+
Step 7
国家高性能计算中心（合肥） 512009-7-24
Systolic乘法
c1,1 c1,2 c1,3 c1,4
c2,1 c2,2 c2,3 c2,4
c3,1 c3,2 c3,3 c3,4
a1,1b1,4
+
a2,1b1,3
+
a3,1b1,2
+
a2,2b2,4
+
a3,2b2,3
+
a3,3b3,4
+
Step 8
国家高性能计算中心（合肥） 522009-7-24
Systolic乘法
c1,1 c1,2 c1,3 c1,4
c2,1 c2,2 c2,3 c2,4
c3,1 c3,2 c3,3 c3,4
a2,1b1,4
+
a3,1b1,3
+
a3,2b2,4
+
Step 9
国家高性能计算中心（合肥） 532009-7-24
Systolic乘法
c1,1 c1,2 c1,3 c1,4
c2,1 c2,2 c2,3 c2,4
c3,1 c3,2 c3,3 c3,4
a3,1b1,4
+
Step 10
国家高性能计算中心（合肥） 542009-7-24
Systolic乘法
P1,1
c1,1
P1,2
c1,2
P1,3
c1,3
P1,4
c1,4
P2,1
c2,1
P2,2
c2,2
P2,3
c2,3
P2,4
c2,4
P3,1
c3,1
P3,2
c3,2
P3,3
c3,3
P3,4
c3,4
Over
c1,1 = a1,1 b1,1 + a1,2 b2,1 + a1,3 b3,1 + a1,4 b4,1
c1,2 = a1,1 b1,2 + a1,2 b2,2 + a1,3 b3,2 + a1,4 b4,2
…………
c3,4 = a3,1 b1,4 + a3,2 b2,4 + a3,3 b3,4 + a3,4 b4,4
国家高性能计算中心（合肥） 552009-7-24
Systolic乘法
Systolic算法
//输入,Am× n,Bn× k; 输出,Cm× k
Begin
for i=1 to m par- do
for j=1 to k par-do
(i) ci,j = 0
(ii) while Pi,j 收到 a和 b时 do
ci,j = ci,j +ab
if i < m then 发送 b给 Pi+1,j endif
if j < k then 发送 a给 Pi,j+1 endif
endwhile
endfor
endfor
End
9.4 矩阵乘法
9.4.1 简单并行分块乘法
9.4.2 Cannon乘法
9.4.3 Fox乘法
9.4.4 Systolic乘法
9.4.5 DNS乘法国家高性能计算中心（合肥） 572009-7-24
DNS乘法
背景,由 Dekel,Nassimi和 Sahni提出的 SIMD-CC上的矩阵乘法,处理器数目为 n3,运行时间为 O(logn),是一种速度很快的算法。
基本思想,通过一到一和一到多的播送办法，使得处理器 (k,i,j)拥有 ai,k和 bk,j,
进行本地相乘,再沿 k方向进行单点积累求和,结果存储在处理器 (0,i,j)中。
处理器编号,处理器数 p=n3= (2q)3=23q,处理器 Pr位于位臵 (k,i,j),
这里 r=kn2+in+j,(0≤i,j,k≤n -1)。位于 (k,i,j)的处理器 Pr的三个寄存器
Ar,Br,Cr分别表示为 A[k,i,j],B[k,i,j]和 C[k,i,j],初始时均为 0。
算法,初始时 ai,j和 bi,j存储于寄存器 A[0,i,j]和 B[0,i,j];
① 数据复制,A,B同时在 k维复制 (一到一播送 )；
A在 j维复制 (一到多播送 ); B在 i维复制 (一到多播送 );
② 相乘运算,所有处理器的 A,B寄存器两两相乘 ;
③ 求和运算,沿 k方向进行单点积累求和 ;
国家高性能计算中心（合肥） 582009-7-24
0
0
0
0
0
0
3
7
1
-5
0
0
4
8
2
-6
1
－ 5
3
7
1
-5
4
8
2
-6
2
－ 6
4
8
3
7
2
－ 5
3
7
1
-5
3
8
1
-6
2
－ 6
4
8
4
7
2
7
3
-5
1
-5
3
-6
1
-6
2
8
4
8
4
7
14
-15-5
16 32
28
139
1410-18-6
(b)(a) (c)
(e)(d) (f)
图 9.12
示例
C00=1× (-5)+2× 7=9
C01=1× (-6)+2× 8=10
C10=3× (-5)+4× 7=13
C11=3× (-6)+4× 8=14
BAC
B
A
求
87
65
43
21
k j
i
初始加载 (b)A,B沿 k维复制 (c)A沿 j维复制
(d)B沿 i维复制 (e)点积 (f)沿 k维求和
B
B
A
A
000 010
001 011
100 110
101 111
P0
P1
P2
P3
P4
P5
P6
P7
国家高性能计算中心（合肥） 592009-7-24
算法描述,
//令 r(m)表示 r的第 m位取反；
//{p,rm=d}表示 r(0≤r≤ p-1)的集合,这里 r的二
//进制第 m位为 d;
//输入,An× n,Bn× n; 输出,Cn× n
Begin //以 n=2,p=8=23举例,q=1,r=(r2r1r0)2
(1)for m=3q-1 to 2q do //按 k维复制 A,B,m=2
for all r in {p,rm=0} par-do //r2=0的 r
(1.1) Ar(m)? Ar //A(100)?A(000)等
(1.2) Br(m)? Br //B(100)?B(000)等
endfor
endfor
(2)for m=q-1 to 0 do //按 j维复制 A,m=0
for all r in {p,rm= r2q+m} par-do //r0=r2的 r
Ar(m)? Ar //A(001)?A(000),A(100)?A(101)
endfor //A(011)?A(010),A(110)?A(111)
endfor
(3)for m=2q-1 to q do //按 i维复制 B,m=1
for all r in {p,rm= rq+m} par-do//r1=r2的 r
Br(m)? Br //B(010)?B(000),B(100)?B(110)
endfor //B(011)?B(001),B(101)?B(111)
endfor
(4)for r=0 to p-1 par-do //相乘,all Pr
Cr=Ar× Br
endfor
(5)for m=2q to 3q-1 do //求和,m=2
for r=0 to p-1 par-do
Cr=Cr+Cr(m)
endfor
endfor
End
DNS乘法国家高性能计算中心（合肥） 602009-7-24
k = 3
0,3 1,3 2,3 3,3
0,2 1,2 2,2 3,2
0,1 1,1 2,1 3,1
1,3 2,3 3,3
0,2 1,2 2,2 3,2
0,1 1,1 2,1 3,1
0,0 1,0 2,0 3,0 0,0 1,0 2,0 3,0
k = 2
k = 1
k = 0
A,B A
0,3
i
j
(b)(a)
k
示例
DNS乘法国家高性能计算中心（合肥） 612009-7-24
k j
i
C [ 0,0 ]
=
A [ 0,3 ] × B [ 3,0 ]
3,3 3,3
3,2
3,1
3,0
0,2 1,2 2,2 3,2
0,1 1,1 2,1 3,1
2,3
2,2
2,1
2,0
1,3
1,2
1,1
1,0
1,0 2,0 3,0
0,0 1,0 2,0 3,0
0,0 1,0 2,0 3,0
0,0 1,0 2,0 3,0
0,3
0,2
0,1
0,0
A
0,3 1,3 2,3 3,3
0,0
( d )( c )
B
0,3 1,3 2,3 3,3
0,3 1,3 2,3 3,3
0,3 1,3 2,3 3,3
0,2 1,2 2,2 3,2
0,2 1,2 2,2 3,2
0,2 1,2 2,2 3,2
0,1 1,1 2,1 3,1
0,1 1,1 2,1 3,1
0,1 1,1 2,1 3,1
+
A [ 0,2 ] × B [ 2,0 ]
+
A [ 0,1 ] × B [ 1,0 ]
+
A [ 0,0 ] × B [ 0,0 ]
3,3
3,3
3,2 3,2 3,2
3,1 3,1 3,1
3,0 3,0 3,0
2,3 2,3 2,3
2,22,2
2,2
2,1 2,1 2,1
2,0 2,0 2,0
1,1
1,1
1,1
1,3 1,3 1,3
1,2
1,2 1,2
1,0 1,0 1,0
0,0
0,0 0,0
0,3 0,3 0,3
0,2 0,2 0,2
0,1 0,1 0,1
DNS乘法
2003年 9月第三篇 并行数值算法第八章 基本通讯操作第九章 稠密矩阵运算第十章 线性方程组的求解第十一章 快速傅里叶变换第九章 稠密矩阵运算
9.1 矩阵的划分
9.2 矩阵转臵
9.3 矩阵 -向量乘法
9.4 矩阵乘法
9.1 矩阵的划分
9.1.1 带状划分
9.1.2 棋盘划分国家高性能计算中心（合肥） 52009-7-24
带状划分
16× 16阶矩阵,p=4
列块带状划分 行循环带状划分
P P P P
0
4
8
1 2
1
5
9
1 3
2
6
1 0
1 4
3
7
1 1
1 5
P
P
P
P
( b )
( a )
图 9,1
0 1 2 3
3
2
1
0
0 1 2 3 4 5 6 7 8 9 1 0 1 1 1 2 1 3 1 4 1 5
国家高性能计算中心（合肥） 62009-7-24
带状划分
示例,p＝ 3,27× 27矩阵的 3种带状划分
9.1 矩阵的划分
9.1.1 带状划分
9.1.2 棋盘划分国家高性能计算中心（合肥） 82009-7-24
棋盘划分
8× 8阶矩阵,p=16
块棋盘划分 循环棋盘划分
( b )( a )
图 9,2
6
1 0
1 4
7
1 1
1 5
4
8
1 2
5
9
1 3
6
1 0
1 4
( 0,0 )
( 1,0 )
( 2,0 )
( 3,0 )
( 5,0 )
( 4,0 )
( 7,0 )
( 6,0 )
( 0,1 )
( 1,1 )
( 2,1 )
( 3,1 )
( 5,1 )
( 4,1 )
( 7,1 )
( 6,1 )
( 0,2 )
( 1,2 )
( 2,2 )
( 3,2 )
( 5,2 )
( 4,2 )
( 7,2 )
( 6,2 )
( 0,3 )
( 1,3 )
( 0,4 )
( 1,4 )
( 2,3 )
( 3,3 )
( 2,4 )
( 3,4 )
( 5,3 )
( 4,3 )
( 5,4 )
( 4,4 )
( 7,3 )
( 6,3 )
( 7,4 )
( 6,4 )
( 0,5 )
( 1,5 )
( 0,6 )
( 1,6 )
( 2,5 )
( 3,5 )
( 2,6 )
( 3,6 )
( 5,5 )
( 4,5 )
( 5,6 )
( 4,6 )
( 7,5 )
( 6,5 )
( 7,6 )
( 6,6 )
( 0,7 )
( 1,7 )
( 2,7 )
( 3,7 )
( 5,7 )
( 4,7 )
( 7,7 )
( 6,7 )
( 0,0 )
( 2,0 )
( 3,0 )
( 5,0 )
( 4,0 )
( 7,0 )
( 6,0 )
( 1,0 )
( 6,1 )
( 1,1 )
( 0,4 )
( 1,4 )
( 2,1 )
( 3,1 )
( 2,4 )
( 3,4 )
( 5,1 )
( 4,1 )
( 5,4 )
( 4,4 )
( 7,1 )( 7,4 )
( 6,4 )
( 0,1 )
( 1,2 )
( 0,5 )
( 1,5 )
( 2,2 )
( 3,2 )
( 2,5 )
( 3,5 )
( 5,2 )
( 4,2 )
( 5,5 )
( 4,5 )
( 7,2 )
( 6,2 )
( 7,5 )
( 6,5 )
( 0,2 )
( 1,3 )
( 0,6 )
( 1,6 )
( 2,3 )( 2,6 )
( 3,6 )
( 5,3 )
( 4,3 )
( 5,6 )
( 4,6 )
( 7,3 )
( 6,3 )
( 7,6 )
( 6,6 )
( 0,3 )
( 3,3 )
( 0,7 )
( 1,7 )
( 2,7 )
( 3,7 )
( 5,7 )
( 4,7 )
( 7,7 )
( 6,7 )
3
7
1 1
1 5
0
4
8
1 2
1
5
9
1 3
P P P P P P P P
2 3 0 1 2
P P P P P P P P
P P P P P P P P
P P P P P P P P
国家高性能计算中心（合肥） 92009-7-24
棋盘划分
示例,p＝ 4,16× 16矩阵的 3种棋盘划分第九章 稠密矩阵运算
9.1 矩阵的划分
9.2 矩阵转臵
9.3 矩阵 -向量乘法
9.4 矩阵乘法
9.2 矩阵转臵
9.2.1 棋盘划分的矩阵转臵
9.2.2 带状划分的矩阵转臵国家高性能计算中心（合肥） 122009-7-24
棋盘划分的矩阵转置
网孔连接
情形 1,p=n2。
通讯步 转臵后
P
( 3,0 )
( 1,0 ) ( 1,2 ) ( 1,3 )
( 2,0 ) ( 2,1 ) ( 2,3 )
P
( 3,1 )
P
( 3,2 )
P P P P
( b )( a )
图 9,3
( 0,0 )
( 0,1 )
( 0,2 )
( 0,3 )
( 1,0 )
( 1,1 )
( 1,2 )
( 1,3 )
( 2,0 )
( 2,1 )
( 2,2 )
( 2,3 )
( 3,0 )
( 3,1 )
( 3,2 )
( 3,3 )
( 0,1 ) ( 0,2 ) ( 0,3 )
P
( 0,0 )
( 1,1 )
( 2,2 )
( 3,3 )
3
7
1 1
1 5
2
6
1 0
1 4
1
5
9
1 31 2
0
4
8
3
7
1 1
1 5
2
6
1 4
1 0
1
9
1 3
5
P P P P P P PP
P P P P P P PP
P P P P P P PP
P
1 2
4
8
0
国家高性能计算中心（合肥） 132009-7-24
棋盘划分的矩阵转置
情形 2,p<n2。
- 划分,
- 算法,① 按 mesh连接进行块转臵 (不同处理器间 )
② 进行块内转臵 (同一处理器内 )
通讯步 转臵后
P P P
P P P
P P P
P P P
( b )
( a )
图 9,4
P P P P
P P P P
P P P P
P P P P
0
4
8
1 2
1
5
9
1 3
2
6
1 0
1 4
3
7
1 1
1 5
( 0,0 )
( 1,0 )
( 0,2 )
( 1,2 )
( 1,4 )
( 0,4 )
( 1,6 )
( 0,6 )
( 0,1 )
( 1,1 )
( 0,3 )
( 1,3 )
( 1,5 )
( 0,5 )
( 1,7 )
( 0,7 )
( 2,0 )
( 3,0 )
( 2,2 )
( 3,2 )
( 3,4 )
( 2,4 )
( 3,6 )
( 2,6 )
( 2,1 )
( 3,1 )
( 2,3 )
( 3,3 )
( 3,5 )
( 2,5 )
( 3,7 )
( 2,7 )
( 4,0 )
( 5,0 )
( 4,2 )
( 5,2 )
( 5,4 )
( 4,4 )
( 5,6 )
( 4,6 )
( 4,1 )
( 5,1 )
( 4,3 )
( 5,3 )
( 5,5 )
( 4,5 )
( 5,7 )
( 4,7 )
( 6,0 )
( 7,0 )
( 6,2 )
( 7,2 )
( 7,4 )
( 6,4 )
( 7,6 )
( 6,6 )
( 6,1 )
( 7,1 )
( 6,3 )
( 7,3 )
( 7,5 )
( 6,5 )
( 7,7 )
( 6,7 )
4
8
1 2
0
P
P
1
5
9
1 3
P
2
6
1 0
1 4
P
3
7
1 1
1 5
( 1,0 )
( 2,0 )
( 3,0 )
( 5,0 )
( 4,0 )
( 7,0 )
( 6,0 )
( 0,0 ) ( 0,1 )
( 1,1 )
( 2,1 )
( 3,1 )
( 5,1 )
( 4,1 )
( 7,1 )
( 6,1 )
( 0,2 )
( 1,2 )
( 2,2 )
( 3,2 )
( 5,2 )
( 4,2 )
( 7,2 )
( 6,2 )
( 0,3 )
( 1,3 )
( 2,3 )
( 3,3 )
( 5,3 )
( 4,3 )
( 7,3 )
( 6,3 )
( 0,4 )
( 1,4 )
( 2,4 )
( 3,4 )
( 5,4 )
( 4,4 )
( 7,4 )
( 6,4 )
( 0,5 )
( 1,5 )
( 2,5 )
( 3,5 )
( 5,5 )
( 4,5 )
( 7,5 )
( 6,5 )
( 0,6 )
( 1,6 )
( 2,6 )
( 3,6 )
( 5,6 )
( 4,6 )
( 7,6 )
( 6,6 )
( 0,7 )
( 1,7 )
( 2,7 )
( 3,7 )
( 5,7 )
( 4,7 )
( 7,7 )
( 6,7 )
子 块划分成p 个大小为A pnpnnn
通讯?)/(2// 2 pnttp ws?
计算pn2// 2
运行时间pntptpnT wsp /222 22
国家高性能计算中心（合肥） 142009-7-24
棋盘划分的矩阵转置
超立方连接
划分,
算法,
①
② 对 Aij递归应用 ① 进行转臵，直至分块矩阵的元素处于同一处理器；
③ 进行同一处理器的内部转臵。
运行时间,
子 块划分成p 个大小为A pnpnnn
p
p
ntt
p
n
p
p
ntt
p
np
p
ntt
p
nT
ws
wswsp
l o g)
2
l o g)2
2
//l o g)2
2
22
2222
（
递归步：，（选路：，内部转置（
2212 21112221 1211 AA AAAA AAA 转置为＝将国家高性能计算中心（合肥） 152009-7-24
棋盘划分的矩阵转置
超立方连接,示例
( 0,5 )
( 1,5 )
P P P
( 2,5 )
( 3,5 )
P
P P P
( 5,5 )
( 4,5 )
P
P P P
( 7,5 )
( 6,5 )
P
( 0,0 ) ( 0,1 )
( 1,0 ) ( 1,1 )
( 0,2 ) ( 0,3 )
( 1,2 ) ( 1,3 )
P
( 0,4 ) ( 0,5 )
( 1,4 ) ( 1,5 )
P
( 0,6 ) ( 0,7 )
( 1,6 ) ( 1,7 )
P
( 2,0 ) ( 2,1 )
( 3,0 ) ( 3,1 )
P
( 2,2 ) ( 2,3 )
( 3,2 ) ( 3,3 )
P
( 2,4 ) ( 2,5 )
( 3,4 ) ( 3,5 )
P
( 2,6 ) ( 2,7 )
( 3,6 ) ( 3,7 )
( 5,0 ) ( 5,1 )
( 4,0 ) ( 4,1 )
( 5,2 ) ( 5,3 )
( 4,2 ) ( 4,3 )
P
( 5,4 ) ( 5,5 )
( 4,4 ) ( 4,5 )
P
( 5,6 ) ( 5,7 )
( 4,6 ) ( 4,7 )
P
( 7,0 ) ( 7,1 )
( 6,0 ) ( 6,1 )
P
( 7,2 ) ( 7,3 )
( 6,2 ) ( 6,3 )
P
( 7,4 ) ( 7,5 )
( 6,4 ) ( 6,5 )
P
( 7,6 ) ( 7,7 )
( 6,6 ) ( 6,7 )
( b )( a )
图 9,6
0
4
8
1 2
1
5
9
2
6
1 0
1 4
3
7
1 1
1 5
0
8
4
1 2
1
9
5
1 3
6
2
1 0
1 4
7
3
1 1
1 5
( 0,0 )
( 1,0 )
( 2,0 )
( 3,0 )
( 5,0 )
( 4,0 )
( 7,0 )
( 6,0 )
( 0,1 )
( 1,1 )
( 2,1 )
( 3,1 )
( 5,1 )
( 4,1 )
( 7,1 )
( 6,1 )
P P P P P P P P
( 0,2 )
( 1,2 )
( 2,2 )
( 3,2 )
( 5,2 )
( 4,2 )
( 7,2 )
( 6,2 )
( 0,3 )
( 1,3 )
( 2,3 )
( 3,3 )
( 5,3 )
( 4,3 )
1 3
( 7,3 )
( 6,3 )
( 0,4 )
( 1,4 )
( 2,4 )
( 3,4 )
( 5,4 )
( 4,4 )
( 7,4 )
( 6,4 )
( 0,6 )
( 1,6 )
( 2,6 )
( 3,6 )
( 5,6 )
( 4,6 )
( 7,6 )
( 6,6 )
( 0,7 )
( 1,7 )
( 2,7 )
( 3,7 )
( 5,7 )
( 4,7 )
( 7,7 )
( 6,7 )
9.2 矩阵转臵
9.2.1 棋盘划分的矩阵转臵
9.2.2 带状划分的矩阵转臵国家高性能计算中心（合肥） 172009-7-24
带状划分的矩阵转置
划分,An× n分成 p个 (n/p)× n大小的带
算法,
① Pi有 p-1个 (n/p)× (n/p)大小子块发送到另外 p-1个处理器中 ;
② 每个处理器本地交换相应的元素
P
P
P
P
n
图 9,7
0
1
2
3
第九章 稠密矩阵运算
9.1 矩阵的划分
9.2 矩阵转臵
9.3 矩阵 -向量乘法
9.4 矩阵乘法
9.3 矩阵 -向量乘法
9.3.1 带状划分的矩阵 -向量乘法
9.3.2 棋盘划分的矩阵 -向量乘法国家高性能计算中心（合肥） 202009-7-24
带状划分的矩阵 -向量乘法
划分 (行带状划分 ),Pi存放 xi和 ai,0,ai,1,…,ai,n-1,并输出 yi
算法,对 p=n情形
① 每个 Pi向其他处理器播送 xi(多到多播送 )；
② 每个 Pi计算；
注,对 p<n情形,算法中 Pi要播送 X中相应的 n/p个分量
(1)超立方连接的计算时间
(2)网孔连接的计算时间充分大时项是多到多的播送时间后项是乘法时间，前
pntpt
p
n
pt
p
npt
p
nT
ws
wsp
//l o g
21//)1(l o g
2
2
充分大时项是多到多的播送时间后项是乘法时间，前
pntpt
p
n
pt
p
ntp
p
nT
ws
wsp
//)1(2
21//)1()1(2
2
2
国家高性能计算中心（合肥） 212009-7-24
带状划分的矩阵 -向量乘法
示例矩 阵 A
0
1
p - 1
n / p
向 量 x 处 理 器
0
1
p - 1
n
( b )( a )
0
0
0
0
0
1
1
1
1
1 p - 1
矩 阵 A
0
1
p - 1
向 量 y
( d )( c )
图 9,8
P
0
P
1
P
p - 1
P
0
P
1
P
p - 1
P
0
P
1
P
p - 1
P
0
P
1
P
p - 1
p - 1
p - 1
p - 1
p - 1
9.3 矩阵 -向量乘法
9.3.1 带状划分的矩阵 -向量乘法
9.3.2 棋盘划分的矩阵 -向量乘法国家高性能计算中心（合肥） 232009-7-24
棋盘划分的矩阵 -向量乘法
划分 (块棋盘划分 ),Pij存放 ai,j,xi臵入 Pi,i中
算法,对 p=n2情形
① 每个 Pi,i向 Pj,i播送 xi(一到多播送 )；
② 按行方向进行乘 -加与积累运算，最后一列 Pi,n-1收集的结果为 yi；
注,对 p<n2情形,p个处理器排成 的二维网孔，
算法中 Pi,i向 Pj,i播送 X中相应的 个分量
(1)网孔连接的计算时间 Tp(CT),
.X中相应分量臵入 Pi,i的通讯时间,
.按列一到多播送时间,
.按行单点积累的时间,ptpt
p
npt
p
nT
hwsp 3loglog
2
pp?
pn/
pttpnt hws
)1(lo g)( ptptpnt hws
)1(lo g)( ptptpnt hws
国家高性能计算中心（合肥） 242009-7-24
棋盘划分的矩阵 -向量乘法
示例矩 阵 A 向 量 x
n
( d )( c )
图 9,9
P
0
P
1
P
n /
矩 阵 A
向 量
y
P
1
P
0
√ p - 1
P
( b )( a )
√ pn /
√ p - 1
√ p
P
p - 1
P
p - 1
P
P
√ p
2 √ p
P
P
√ p
2 √ p
国家高性能计算中心（合肥） 252009-7-24
带状与棋盘划分比较
以网孔链接为例
网孔上带状划分的运行时间
网孔上棋盘划分的运行时间
棋盘划分要比带状划分快。
)5.9()1(2
2
wsp ntptp
nT
)6.9(3loglog
2
ptptpnptpnT hwsp
第九章 稠密矩阵运算
9.1 矩阵的划分
9.2 矩阵转臵
9.3 矩阵 -向量乘法
9.4 矩阵乘法
9.4 矩阵乘法
9.4.1 简单并行分块乘法
9.4.2 Cannon乘法
9.4.3 Fox乘法
9.4.4 Systolic乘法
9.4.5 DNS乘法国家高性能计算中心（合肥） 282009-7-24
矩阵乘法符号及定义
1,11,10,1
1,11,10,1
1,01,00,0
1,11,10,1
1,11,10,1
1,01,00,0
1,11,10,1
1,11,10,1
1,01,00,0
,)()()(
nnnn
n
n
nnnn
n
n
nnnn
n
n
nnijnnijnnij
bbb
bbb
bbb
aaa
aaa
aaa
ccc
ccc
ccc
BACcCbBaA
设
j
i
A BC
1
0
n
k
kjikij bac
A中元素的第 1下标与 B中元素的第 2下标相一致（对准）
国家高性能计算中心（合肥） 292009-7-24
矩阵乘法并行实现方法
计算结构：二维阵列
空间对准 (元素已加载到阵列中）
Cannon’s,Fox’s,DNS
时间对准 (元素未加载到阵列中）
Systolic A0,0
B0,0
A1,0
B1,0
A2,0
B2,0
A3,0
B3,0
A0,1
B0,1
A1,1
B1,1
A2,1
B2,1
A3,1
B3,1
A0,2
B0,2
A1,2
B1,2
A2,2
B2,2
A3,2
B3,2
A0,3
B0,3
A1,3
B1,3
A2,3
B2,3
A3,3
B3,3
国家高性能计算中心（合肥） 302009-7-24
简单并行分块乘法
分块,A,B和 C分成 的方块阵 Ai,j,Bi,j和 Ci,j,大小均为
p个处理器编号为,Pi,j存放 Ai,j,Bi,j和 Ci,j。
算法,
① 通讯,每行处理器进行 A矩阵块的多到多播送 (得到 Ai,k,k=0~ )
每列处理器进行 B矩阵块的多到多播送 (得到 Bk,j,k=0~ )
② 乘 -加运算,Pi,j做
运行时间
(1)超立方连接：
① 的时间
② 的时间
ppp pnpn?
),.,,,,.,,,( 1,11,00,0 ppp PPP
1－p
10pk kjikij BAC
))1(lo g(2 21 ppntptt ws
pnpt pn /)( 332
1－p
国家高性能计算中心（合肥） 312009-7-24
简单并行分块乘法
运行时间
(1)超立方连接：
(2)二维环绕网孔连接：
① 的时间：
② 的时间 t2=n3/p
注
(1)本算法的缺点是对处理器的存储要求过大每个处理器有 个块，每块大小为 n2/p,
所以需要,p个处理器共需要,
是串行算法的 倍
(2)p=n2时,t(n)=O(n),c(n)=O(n3)
p
ntptpt
p
ntt
wsws
22
1 22)1((2 ）
p
ntpt
p
nT
wsp
23 22
p2
)/( 2 pnO )( 2 pnO
p
9.4 矩阵乘法
9.4.1 简单并行分块乘法
9.4.2 Cannon乘法
9.4.3 Fox乘法
9.4.4 Systolic乘法
9.4.5 DNS乘法国家高性能计算中心（合肥） 332009-7-24
Cannon乘法
分块,A,B和 C分成 的方块阵 Ai,j,Bi,j和 Ci,j,大小均为 p个处理器编号为,Pi,j存放
Ai,j,Bi,j和 Ci,j(n > > p)
ppp
pnpn? ),.,,,,.,,,( 1,11,00,0 ppp PPP
P0,0
P1,0
P2,0
P3,0
P0,1
P1,1
P2,1
P3,1
P0,2
P1,2
P2,2
P3,2
P0,3
P1,3
P2,3
P3,3
p
n
n
p
国家高性能计算中心（合肥） 342009-7-24
Cannon乘法
算法原理 (非形式描述 )
① 所有块 Ai,j(0≤i,j≤ ) 向左循环移动 i步 (按行移位 )；
所有块 Bi,j(0≤i,j≤ ) 向上循环移动 j步 (按列移位 )；
② 所有处理器 Pi,j做执行 Ai,j和 Bi,j的乘 -加运算 ;
③ A的每个块向左循环移动一步 ;
B的每个块向上循环移动一步 ;
④ 转 ② 执行 次 ;
1－p
1－p
1－p
国家高性能计算中心（合肥） 352009-7-24
Cannon乘法
示例,A4× 4,B4× 4,p=16
A0,0
A1,0
A2,0
A3,0
A0,1
A1,1
A2,1
A3,1
A0,2
A1,2
A2,2
A3,2
A0,3
A1,3
A2,3
A3,3
B0,0
B1,0
B2,0
B3,0
B0,1
B1,1
B2,1
B3,1
B0,2
B1,2
B2,2
B3,2
B0,3
B1,3
B2,3
B3,3
Initial alignment of A Initial alignment of B
国家高性能计算中心（合肥） 362009-7-24
Cannon乘法
示例,A4× 4,B4× 4,p=16
A and B after initial alignment and shifts after every step
A0,0
B0,0
A1,1
B1,0
A2,2
B2,0
A3,3
B3,0
A0,1
B1,1
A1,2
B2,1
A2,3
B3,1
A3,0
B0,1
A0,2
B2,2
A1,3
B3,2
A2,0
B0,2
A3,1
B1,2
A0,3
B3,3
A1,0
B0,3
A2,1
B1,3
A3,2
B2,3
国家高性能计算中心（合肥） 372009-7-24
Cannon乘法
示例,A4× 4,B4× 4,p=16
After first shift
A0,1
B1,0
A1,2
B2,0
A2,3
B3,0
A3,0
B0,0
A0,2
B2,1
A1,3
B3,1
A2,0
B0,1
A3,1
B3,1
A0,3
B3,2
A1,0
B0,2
A2,1
B1,2
A3,2
B2,2
A0,0
B0,3
A1,1
B1,3
A2,2
B2,3
A3,3
B3,3
After second shift
A0,2
B2,0
A1,3
B3,0
A2,0
B0,0
A3,1
B1,0
A0,3
B3,1
A1,0
B0,1
A2,1
B1,1
A3,2
B2,1
A0,0
B0,2
A1,1
B1,2
A2,2
B2,2
A3,3
B3,2
A0,1
B1,3
A1,2
B2,3
A2,3
B3,3
A3,0
B0,3
After third shift
A0,3
B3,0
A1,0
B0,0
A2,1
B1,0
A3,2
B2,0
A0,0
B0,1
A1,1
B1,1
A2,2
B2,1
A3,3
B3,1
A0,1
B1,2
A1,2
B2,2
A2,3
B3,2
A3,0
B0,2
A0,2
B2,3
A1,3
B3,3
A2,0
B0,3
A3,1
B1,3
国家高性能计算中心（合肥） 382009-7-24
Cannon乘法
算法描述,Cannon分块乘法算法
//输入,An× n,Bn× n; 输出,Cn× n
Begin
(1)for k=0 to do
for all Pi,j par-do
(i) if i>k then
Ai,j? Ai,(j+1)mod
endif
(ii)if j>k then
Bi,j? B(i+1)mod,j
endif
endfor
endfor
(2)for all Pi,j par-do Ci,j=0 endfor
(3)for k=0 to do
for all Pi,j par-do
(i) Ci,j=Ci,j+Ai,jBi,j
(ii) Ai,j? Ai,(j+1)mod
(iii)Bi,j? B(i+1)mod,j
endfor
endfor
End
时间分析：
)/(
))/(()1()(
)(
3
3
321
pnO
pnpOOpO
TTTnT p
1－p
p
p
p
p
1－p
9.4 矩阵乘法
9.4.1 简单并行分块乘法
9.4.2 Cannon乘法
9.4.3 Fox乘法
9.4.4 Systolic乘法
9.4.5 DNS乘法国家高性能计算中心（合肥） 402009-7-24
Fox乘法
分块,同 Cannon分块算法
算法原理
① Ai,i向所在行的其他处理器进行一到多播送；
② 各处理器将收到的 A块与原有的 B块进行乘 -加运算 ;
③ B块向上循环移动一步 ;
④ 如果 Ai,j是上次第 i行播送的块，本次选择 向所在行的其他处理器进行一到多播送 ;
⑤ 转 ② 执行 次 ;1?p
pjiA m od)1(,?
A0,0
B0,0
A1,0
B1,0
A2,0
B2,0
A3,0
B3,0
A0,1
B0,1
A1,1
B1,1
A2,1
B2,1
A3,1
B3,1
A0,2
B0,2
A1,2
B1,2
A2,2
B2,2
A3,2
B3,2
A0,3
B0,3
A1,3
B1,3
A2,3
B2,3
A3,3
B3,3
国家高性能计算中心（合肥） 412009-7-24
Fox乘法
示例,A4× 4,B4× 4,p=16
(a) (b)
A0,0
B0,0
B1,0
B2,0
B3,0
B0,1
A1,1
B1,1
B2,1
B3,1
B0,2
B1,2
A2,2
B2,2
B3,2
B0,3
B1,3
B2,3
A3,3
B3,3
B1,0
B2,0
B3,0
A0,1
B1,1
B2,1
B3,1
B0,1
B1,2
B3,2
B0,2
B1,3
B2,3
B0,3
A1,2
B2,2
A2,3
B3,3
A3,0
B0,0
国家高性能计算中心（合肥） 422009-7-24
Fox乘法
示例,A4× 4,B4× 4,p=16
(c) (d)
B2,0
B3,0
B2,1
B3,1
B0,1
B3,2
B0,2
B1,2
B2,3
B0,3
B1,3
B3,0
B1,0
B3,1
B0,1
B2,1
B3,2
B1,2
B0,3
B2,3
B0,2
B1,3
B2,0
A0,2
B2,2
A1,3
B3,3
A2,0
B0,0
B1,0
A3,1
B1,1
A0,3
B3,3
A1,0
B0,2
A2,1
B1,1
A3,2
B2,2
9.4 矩阵乘法
9.4.1 简单并行分块乘法
9.4.2 Cannon乘法
9.4.3 Fox乘法
9.4.4 Systolic乘法
9.4.5 DNS乘法国家高性能计算中心（合肥） 442009-7-24
Systolic乘法
a1,4
b4,1
b3,1
b2,1
b2,2
b4,2
b3,2
b2,3
b3,3
b4,3
b2,4
b3,4
b4,4
a1,3a1,1 a1,2
a2,4a2,1 a2,2 a2,3
a3,1 a3,2 a3,3 a3,4
b1,1
b1,2
b1,3
b1,4
Step 1
P1,1
c1,1
P1,2
c1,2
P1,3
c1,3
P1,4
c1,4
P2,1
c2,1
P2,2
c2,2
P2,3
c2,3
P2,4
c2,4
P3,1
c3,1
P3,2
c3,2
P3,3
c3,3
P3,4
c3,4
国家高性能计算中心（合肥） 452009-7-24
Systolic乘法
c1,1 c1,2 c1,3 c1,4
c2,1 c2,2 c2,3 c2,4
c3,1 c3,2 c3,3 c3,4
b3,1
b2,1
b2,2
b4,2
b3,2
b2,3
b3,3
b4,3
b2,4
b3,4
b4,4
a1,3a1,1 a1,2
a2,4a2,1 a2,2 a2,3
a3,1 a3,2 a3,3 a3,4
b1,1
b1,2
b1,3
b1,4
a1,4b4,1
+
Step 2
国家高性能计算中心（合肥） 462009-7-24
Systolic乘法
c1,1 c1,2 c1,3 c1,4
c2,1 c2,2 c2,3 c2,4
c3,1 c3,2 c3,3 c3,4
b2,1
b2,2
b3,2
b2,3
b3,3
b4,3
b2,4
b3,4
b4,4
a1,1 a1,2
a2,1 a2,2 a2,3
a3,1 a3,2 a3,3 a3,4
b1,1
b1,2
b1,3
b1,4
a1,3b3,1
+
a1,4b4,2
+
a2,4b4,1
+
Step 3
国家高性能计算中心（合肥） 472009-7-24
Systolic乘法
c1,1 c1,2 c1,3 c1,4
c2,1 c2,2 c2,3 c2,4
c3,1 c3,2 c3,3 c3,4
b2,2
b2,3
b3,3
b2,4
b3,4
b4,4
a1,1
a2,1 a2,2
a3,1 a3,2 a3,3
b1,1
b1,2
b1,3
b1,4
a1,2b2,1
+
a1,3b3,2
+
a2,3b3,1
+
a1,4b4,3
+
a3,4b4,1
+
a2,4b4,2
+
Step 4
国家高性能计算中心（合肥） 482009-7-24
Systolic乘法
c1,1 c1,2 c1,3 c1,4
c2,1 c2,2 c2,3 c2,4
c3,1 c3,2 c3,3 c3,4
b2,3
b2,4
b3,4
a2,1
a3,1 a3,2
b1,2
b1,3
b1,4
a1,1b1,1
+
a1,2b2,2
+
a2,2b2,1
+
a1,3b3,3
+
a3,3b3,1
+
a2,3b3,2
+
a1,4b4,4
+
a2,4b4,3
+
a3,4b4,2
+
Step 5
国家高性能计算中心（合肥） 492009-7-24
Systolic乘法
c1,1 c1,2 c1,3 c1,4
c2,1 c2,2 c2,3 c2,4
c3,1 c3,2 c3,3 c3,4
b2,4
a3,1
b1,3
b1,4
a1,1b1,2
+
a2,1b1,1
+
a1,2b2,3
+
a3,2b2,1
+
a2,2b2,2
+
a1,3b3,4
+
a2,3b3,3
+
a3,3b3,2
+
a2,4b4,4
+
a3,4b4,3
+
Step 6
国家高性能计算中心（合肥） 502009-7-24
Systolic乘法
c1,1 c1,2 c1,3 c1,4
c2,1 c2,2 c2,3 c2,4
c3,1 c3,2 c3,3 c3,4
b1,4
a1,1b1,3
+
a3,1b1,1
+
a2,1b1,2
+
a1,2b2,4
+
a2,2b2,3
+
a3,2b3,2
+
a2,3b3,4
+
a3,3b3,3
+
a3,4b4,4
+
Step 7
国家高性能计算中心（合肥） 512009-7-24
Systolic乘法
c1,1 c1,2 c1,3 c1,4
c2,1 c2,2 c2,3 c2,4
c3,1 c3,2 c3,3 c3,4
a1,1b1,4
+
a2,1b1,3
+
a3,1b1,2
+
a2,2b2,4
+
a3,2b2,3
+
a3,3b3,4
+
Step 8
国家高性能计算中心（合肥） 522009-7-24
Systolic乘法
c1,1 c1,2 c1,3 c1,4
c2,1 c2,2 c2,3 c2,4
c3,1 c3,2 c3,3 c3,4
a2,1b1,4
+
a3,1b1,3
+
a3,2b2,4
+
Step 9
国家高性能计算中心（合肥） 532009-7-24
Systolic乘法
c1,1 c1,2 c1,3 c1,4
c2,1 c2,2 c2,3 c2,4
c3,1 c3,2 c3,3 c3,4
a3,1b1,4
+
Step 10
国家高性能计算中心（合肥） 542009-7-24
Systolic乘法
P1,1
c1,1
P1,2
c1,2
P1,3
c1,3
P1,4
c1,4
P2,1
c2,1
P2,2
c2,2
P2,3
c2,3
P2,4
c2,4
P3,1
c3,1
P3,2
c3,2
P3,3
c3,3
P3,4
c3,4
Over
c1,1 = a1,1 b1,1 + a1,2 b2,1 + a1,3 b3,1 + a1,4 b4,1
c1,2 = a1,1 b1,2 + a1,2 b2,2 + a1,3 b3,2 + a1,4 b4,2
…………
c3,4 = a3,1 b1,4 + a3,2 b2,4 + a3,3 b3,4 + a3,4 b4,4
国家高性能计算中心（合肥） 552009-7-24
Systolic乘法
Systolic算法
//输入,Am× n,Bn× k; 输出,Cm× k
Begin
for i=1 to m par- do
for j=1 to k par-do
(i) ci,j = 0
(ii) while Pi,j 收到 a和 b时 do
ci,j = ci,j +ab
if i < m then 发送 b给 Pi+1,j endif
if j < k then 发送 a给 Pi,j+1 endif
endwhile
endfor
endfor
End
9.4 矩阵乘法
9.4.1 简单并行分块乘法
9.4.2 Cannon乘法
9.4.3 Fox乘法
9.4.4 Systolic乘法
9.4.5 DNS乘法国家高性能计算中心（合肥） 572009-7-24
DNS乘法
背景,由 Dekel,Nassimi和 Sahni提出的 SIMD-CC上的矩阵乘法,处理器数目为 n3,运行时间为 O(logn),是一种速度很快的算法。
基本思想,通过一到一和一到多的播送办法，使得处理器 (k,i,j)拥有 ai,k和 bk,j,
进行本地相乘,再沿 k方向进行单点积累求和,结果存储在处理器 (0,i,j)中。
处理器编号,处理器数 p=n3= (2q)3=23q,处理器 Pr位于位臵 (k,i,j),
这里 r=kn2+in+j,(0≤i,j,k≤n -1)。位于 (k,i,j)的处理器 Pr的三个寄存器
Ar,Br,Cr分别表示为 A[k,i,j],B[k,i,j]和 C[k,i,j],初始时均为 0。
算法,初始时 ai,j和 bi,j存储于寄存器 A[0,i,j]和 B[0,i,j];
① 数据复制,A,B同时在 k维复制 (一到一播送 )；
A在 j维复制 (一到多播送 ); B在 i维复制 (一到多播送 );
② 相乘运算,所有处理器的 A,B寄存器两两相乘 ;
③ 求和运算,沿 k方向进行单点积累求和 ;
国家高性能计算中心（合肥） 582009-7-24
0
0
0
0
0
0
3
7
1
-5
0
0
4
8
2
-6
1
－ 5
3
7
1
-5
4
8
2
-6
2
－ 6
4
8
3
7
2
－ 5
3
7
1
-5
3
8
1
-6
2
－ 6
4
8
4
7
2
7
3
-5
1
-5
3
-6
1
-6
2
8
4
8
4
7
14
-15-5
16 32
28
139
1410-18-6
(b)(a) (c)
(e)(d) (f)
图 9.12
示例
C00=1× (-5)+2× 7=9
C01=1× (-6)+2× 8=10
C10=3× (-5)+4× 7=13
C11=3× (-6)+4× 8=14
BAC
B
A
求
87
65
43
21
k j
i
初始加载 (b)A,B沿 k维复制 (c)A沿 j维复制
(d)B沿 i维复制 (e)点积 (f)沿 k维求和
B
B
A
A
000 010
001 011
100 110
101 111
P0
P1
P2
P3
P4
P5
P6
P7
国家高性能计算中心（合肥） 592009-7-24
算法描述,
//令 r(m)表示 r的第 m位取反；
//{p,rm=d}表示 r(0≤r≤ p-1)的集合,这里 r的二
//进制第 m位为 d;
//输入,An× n,Bn× n; 输出,Cn× n
Begin //以 n=2,p=8=23举例,q=1,r=(r2r1r0)2
(1)for m=3q-1 to 2q do //按 k维复制 A,B,m=2
for all r in {p,rm=0} par-do //r2=0的 r
(1.1) Ar(m)? Ar //A(100)?A(000)等
(1.2) Br(m)? Br //B(100)?B(000)等
endfor
endfor
(2)for m=q-1 to 0 do //按 j维复制 A,m=0
for all r in {p,rm= r2q+m} par-do //r0=r2的 r
Ar(m)? Ar //A(001)?A(000),A(100)?A(101)
endfor //A(011)?A(010),A(110)?A(111)
endfor
(3)for m=2q-1 to q do //按 i维复制 B,m=1
for all r in {p,rm= rq+m} par-do//r1=r2的 r
Br(m)? Br //B(010)?B(000),B(100)?B(110)
endfor //B(011)?B(001),B(101)?B(111)
endfor
(4)for r=0 to p-1 par-do //相乘,all Pr
Cr=Ar× Br
endfor
(5)for m=2q to 3q-1 do //求和,m=2
for r=0 to p-1 par-do
Cr=Cr+Cr(m)
endfor
endfor
End
DNS乘法国家高性能计算中心（合肥） 602009-7-24
k = 3
0,3 1,3 2,3 3,3
0,2 1,2 2,2 3,2
0,1 1,1 2,1 3,1
1,3 2,3 3,3
0,2 1,2 2,2 3,2
0,1 1,1 2,1 3,1
0,0 1,0 2,0 3,0 0,0 1,0 2,0 3,0
k = 2
k = 1
k = 0
A,B A
0,3
i
j
(b)(a)
k
示例
DNS乘法国家高性能计算中心（合肥） 612009-7-24
k j
i
C [ 0,0 ]
=
A [ 0,3 ] × B [ 3,0 ]
3,3 3,3
3,2
3,1
3,0
0,2 1,2 2,2 3,2
0,1 1,1 2,1 3,1
2,3
2,2
2,1
2,0
1,3
1,2
1,1
1,0
1,0 2,0 3,0
0,0 1,0 2,0 3,0
0,0 1,0 2,0 3,0
0,0 1,0 2,0 3,0
0,3
0,2
0,1
0,0
A
0,3 1,3 2,3 3,3
0,0
( d )( c )
B
0,3 1,3 2,3 3,3
0,3 1,3 2,3 3,3
0,3 1,3 2,3 3,3
0,2 1,2 2,2 3,2
0,2 1,2 2,2 3,2
0,2 1,2 2,2 3,2
0,1 1,1 2,1 3,1
0,1 1,1 2,1 3,1
0,1 1,1 2,1 3,1
+
A [ 0,2 ] × B [ 2,0 ]
+
A [ 0,1 ] × B [ 1,0 ]
+
A [ 0,0 ] × B [ 0,0 ]
3,3
3,3
3,2 3,2 3,2
3,1 3,1 3,1
3,0 3,0 3,0
2,3 2,3 2,3
2,22,2
2,2
2,1 2,1 2,1
2,0 2,0 2,0
1,1
1,1
1,1
1,3 1,3 1,3
1,2
1,2 1,2
1,0 1,0 1,0
0,0
0,0 0,0
0,3 0,3 0,3
0,2 0,2 0,2
0,1 0,1 0,1
DNS乘法
