中国科学技术大学计算机科学与技术系 国家高性能计算中心(合肥) 2004年12月 并 行 计 算 中国科学技术大学计算机科学与技术系 国家高性能计算中心(合肥) 2004年12月

第三篇 并行数值算法 第八章 基本通讯操作 第九章 稠密矩阵运算 第十章 线性方程组的求解 第十一章 快速傅里叶变换 第三篇 并行数值算法 第八章 基本通讯操作 第九章 稠密矩阵运算 第十章 线性方程组的求解 第十一章 快速傅里叶变换

第九章 稠密矩阵运算 9.1 矩阵的划分 9.2 矩阵转置 9.3 矩阵-向量乘法 9.4 矩阵乘法 第九章 稠密矩阵运算 9.1 矩阵的划分 9.2 矩阵转置 9.3 矩阵-向量乘法 9.4 矩阵乘法

9.1 矩阵的划分 9.1.1 带状划分 9.1.2 棋盘划分

带状划分 16×16阶矩阵，p=4 列块带状划分 行循环带状划分 国家高性能计算中心（合肥） 2019/1/3

带状划分 示例：p＝3，27× 27矩阵的3种带状划分 国家高性能计算中心（合肥） 2019/1/3

9.1 矩阵的划分 9.1.1 带状划分 9.1.2 棋盘划分

棋盘划分 8×8阶矩阵，p=16 块棋盘划分 循环棋盘划分 国家高性能计算中心（合肥） 2019/1/3

棋盘划分 示例：p＝4，16×16矩阵的3种棋盘划分 国家高性能计算中心（合肥） 2019/1/3

第九章 稠密矩阵运算 9.1 矩阵的划分 9.2 矩阵转置 9.3 矩阵-向量乘法 9.4 矩阵乘法 第九章 稠密矩阵运算 9.1 矩阵的划分 9.2 矩阵转置 9.3 矩阵-向量乘法 9.4 矩阵乘法

9.2 矩阵转置 9.2.1 棋盘划分的矩阵转置 9.2.2 带状划分的矩阵转置 9.2 矩阵转置 9.2.1 棋盘划分的矩阵转置 9.2.2 带状划分的矩阵转置

棋盘划分的矩阵转置 网孔连接 情形1: p=n2。 通讯步 转置后 国家高性能计算中心（合肥） 2019/1/3

棋盘划分的矩阵转置 情形2: p<n2。 - 算法: ①按mesh连接进行块转置(不同处理器间) ②进行块内转置(同一处理器内) - 划分: - 算法: ①按mesh连接进行块转置(不同处理器间) ②进行块内转置(同一处理器内) 通讯步 转置后 国家高性能计算中心（合肥） 2019/1/3

棋盘划分的矩阵转置 超立方连接 划分: 算法: ① 运行时间: ②对Aij递归应用①进行转置，直至分块矩阵的元素处于同一处理器； ③进行同一处理器的内部转置。 运行时间: 国家高性能计算中心（合肥） 2019/1/3

棋盘划分的矩阵转置 超立方连接：示例 国家高性能计算中心（合肥） 2019/1/3

9.2 矩阵转置 9.2.1 棋盘划分的矩阵转置 9.2.2 带状划分的矩阵转置 9.2 矩阵转置 9.2.1 棋盘划分的矩阵转置 9.2.2 带状划分的矩阵转置

带状划分的矩阵转置 划分: An×n分成p个(n/p)×n大小的带 算法: ①Pi有p-1个(n/p)×(n/p)大小子块发送到另外p-1个处理器中; ②每个处理器本地交换相应的元素 国家高性能计算中心（合肥） 2019/1/3

第九章 稠密矩阵运算 9.1 矩阵的划分 9.2 矩阵转置 9.3 矩阵-向量乘法 9.4 矩阵乘法 第九章 稠密矩阵运算 9.1 矩阵的划分 9.2 矩阵转置 9.3 矩阵-向量乘法 9.4 矩阵乘法

9.3 矩阵-向量乘法 9.3.1 带状划分的矩阵-向量乘法 9.3.2 棋盘划分的矩阵-向量乘法 9.3 矩阵-向量乘法 9.3.1 带状划分的矩阵-向量乘法 9.3.2 棋盘划分的矩阵-向量乘法

带状划分的矩阵-向量乘法 划分(行带状划分): Pi存放xi和ai,0,ai,1,…,ai,n-1, 并输出yi 算法: 对p=n情形 注: 对p<n情形,算法中Pi要播送X中相应的n/p个分量 (1)超立方连接的计算时间 (2)网孔连接的计算时间 国家高性能计算中心（合肥） 2019/1/3

带状划分的矩阵-向量乘法 示例 国家高性能计算中心（合肥） 2019/1/3

9.3 矩阵-向量乘法 9.3.1 带状划分的矩阵-向量乘法 9.3.2 棋盘划分的矩阵-向量乘法 9.3 矩阵-向量乘法 9.3.1 带状划分的矩阵-向量乘法 9.3.2 棋盘划分的矩阵-向量乘法

棋盘划分的矩阵-向量乘法 划分(块棋盘划分): 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的通讯时间: .按列一到多播送时间: .按行单点积累的时间: 国家高性能计算中心（合肥） 2019/1/3

棋盘划分的矩阵-向量乘法 示例 国家高性能计算中心（合肥） 2019/1/3

带状与棋盘划分比较 以网孔链接为例 棋盘划分要比带状划分快。 网孔上带状划分的运行时间 网孔上棋盘划分的运行时间 国家高性能计算中心（合肥） 2019/1/3

第九章 稠密矩阵运算 9.1 矩阵的划分 9.2 矩阵转置 9.3 矩阵-向量乘法 9.4 矩阵乘法 第九章 稠密矩阵运算 9.1 矩阵的划分 9.2 矩阵转置 9.3 矩阵-向量乘法 9.4 矩阵乘法

9. 4 矩阵乘法 9. 4. 1 简单并行分块乘法 9. 4. 2 Cannon乘法 9. 4. 3 Fox乘法 9. 4 9.4 矩阵乘法 9.4.1 简单并行分块乘法 9.4.2 Cannon乘法 9.4.3 Fox乘法 9.4.4 Systolic乘法 9.4.5 DNS乘法

矩阵乘法符号及定义 j i A B C A中元素的第1下标与B中元素的第2下标相一致（对准） 国家高性能计算中心（合肥） 2019/1/3

矩阵乘法并行实现方法 计算结构：二维阵列 空间对准(元素已加载到阵列中） 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 国家高性能计算中心（合肥） 2019/1/3

简单并行分块乘法 分块: 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)超立方连接： ①的时间 ②的时间 国家高性能计算中心（合肥） 2019/1/3

简单并行分块乘法 运行时间 注 (1)超立方连接： (2)二维环绕网孔连接： ①的时间： ②的时间t2=n3/p (1)本算法的缺点是对处理器的存储要求过大 每个处理器有 个块，每块大小为n2/p, 所以需要 , p个处理器共需要 ， 是串行算法的 倍 (2)p=n2时，t(n)=O(n), c(n)=O(n3) 国家高性能计算中心（合肥） 2019/1/3

9. 4 矩阵乘法 9. 4. 1 简单并行分块乘法 9. 4. 2 Cannon乘法 9. 4. 3 Fox乘法 9. 4 9.4 矩阵乘法 9.4.1 简单并行分块乘法 9.4.2 Cannon乘法 9.4.3 Fox乘法 9.4.4 Systolic乘法 9.4.5 DNS乘法

Cannon乘法 分块: A、B和C分成 的方块阵Ai,j、Bi,j和Ci,j, 大小 均为p个处理器编号为 , Pi,j存放Ai,j、Bi,j和Ci,j(n > > p) 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 国家高性能计算中心（合肥） 2019/1/3

Cannon乘法 算法原理 (非形式描述) 所有块Bi,j(0≤i,j≤ )向上循环移动j步(按列移位)； ①所有块Ai,j(0≤i,j≤ )向左循环移动i步(按行移位)； 所有块Bi,j(0≤i,j≤ )向上循环移动j步(按列移位)； ②所有处理器Pi,j做执行Ai,j和Bi,j的乘-加运算; ③A的每个块向左循环移动一步; B的每个块向上循环移动一步; ④转②执行 次; 国家高性能计算中心（合肥） 2019/1/3

Cannon乘法 示例： A4×4， B4×4, p=16 Initial alignment of A Initial alignment of B 国家高性能计算中心（合肥） 2019/1/3

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 国家高性能计算中心（合肥） 2019/1/3

Cannon乘法 示例： A4×4， B4×4, p=16 After first shift After second shift After third shift 国家高性能计算中心（合肥） 2019/1/3

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 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 End 时间分析： 国家高性能计算中心（合肥） 2019/1/3

9. 4 矩阵乘法 9. 4. 1 简单并行分块乘法 9. 4. 2 Cannon乘法 9. 4. 3 Fox乘法 9. 4 9.4 矩阵乘法 9.4.1 简单并行分块乘法 9.4.2 Cannon乘法 9.4.3 Fox乘法 9.4.4 Systolic乘法 9.4.5 DNS乘法

Fox乘法 分块: 同Cannon分块算法 算法原理 ①Ai,i向所在行的其他处理器 进行一到多播送； ②各处理器将收到的A块与原 有的B块进行乘-加运算; ③B块向上循环移动一步; ④如果Ai,j是上次第i行播送的块，本次选择 向所 在行的其他处理器进行一到多播送; ⑤转②执行 次; 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 国家高性能计算中心（合肥） 2019/1/3

Fox乘法 示例： A4×4， B4×4, p=16 (a) (b) A0,0 B0,0 B1,0 B2,0 B3,0 B0,1 A1,1 国家高性能计算中心（合肥） 2019/1/3

Fox乘法 示例： A4×4， B4×4, p=16 (c) (d) B2,0 B2,1 A0,2 B2,2 B2,3 B3,0 B3,1 国家高性能计算中心（合肥） 2019/1/3

9. 4 矩阵乘法 9. 4. 1 简单并行分块乘法 9. 4. 2 Cannon乘法 9. 4. 3 Fox乘法 9. 4 9.4 矩阵乘法 9.4.1 简单并行分块乘法 9.4.2 Cannon乘法 9.4.3 Fox乘法 9.4.4 Systolic乘法 9.4.5 DNS乘法

Systolic乘法 Step 1 b1,4 b1,3 b2,4 b1,2 b2,3 b3,4 b1,1 b2,2 b3,3 b4,4 a1,1 a1,2 a1,3 a1,4 P2, 1 c2,1 P2,2 c2,2 P2,3 c2,3 P2,4 c2,4 a2,1 a2,2 a2,3 a2,4 P3,1 c3,1 P3,2 c3,2 P3,3 c3,3 P3,4 c3,4 a3,1 a3,2 a3,3 a3,4 国家高性能计算中心（合肥） 2019/1/3

Systolic乘法 Step 2 b1,4 b1,3 b2,4 b1,2 b2,3 b3,4 b1,1 b2,2 b3,3 b4,4 a1,4 b4,1 c1,2 c1,3 c1,4 a1,1 a1,2 a1,3 + c2,1 c2,2 c2,3 c2,4 a2,1 a2,2 a2,3 a2,4 c3,1 c3,2 c3,3 c3,4 a3,1 a3,2 a3,3 a3,4 国家高性能计算中心（合肥） 2019/1/3

Systolic乘法 Step 3 b1,4 b1,3 b2,4 b1,2 b2,3 b3,4 b1,1 b2,2 b3,3 b4,4 a1,3 b3,1 c1,2 a1,4 b4,2 c1,3 c1,4 a1,1 a1,2 + + c2,1 a2,4 b4,1 c2,2 c2,3 c2,4 a2,1 a2,2 a2,3 + c3,1 c3,2 c3,3 c3,4 a3,1 a3,2 a3,3 a3,4 国家高性能计算中心（合肥） 2019/1/3

Systolic乘法 Step 4 b1,4 b1,3 b2,4 b1,2 b2,3 b3,4 b1,1 b2,2 b3,3 b4,4 a1,2 b2,1 a1,3 c1,2 b3,2 c1,3 a1,4 b4,3 c1,4 a1,1 + + + c2,1 a2,3 b3,1 c2,2 a2,4 b4,2 c2,3 c2,4 a2,1 a2,2 + + c3,1 a3,4 b4,1 c3,2 c3,3 c3,4 a3,1 a3,2 a3,3 + 国家高性能计算中心（合肥） 2019/1/3

Systolic乘法 Step 5 b1,4 b1,3 b2,4 b1,2 b2,3 b3,4 c1,1 a1,1 b1,1 a1,2 + + + + c2,1 a2,2 b2,1 c2,2 a2,3 b3,2 c2,3 a2,4 b4,3 c2,4 a2,1 + + + c3,1 a3,3 b3,1 c3,2 a3,4 b4,2 c3,3 c3,4 a3,1 a3,2 + + 国家高性能计算中心（合肥） 2019/1/3

Systolic乘法 Step 6 b1,4 b1,3 b2,4 c1,1 c1,2 a1,1 b1,2 c1,3 a1,2 b2,3 + + + c2,1 a2,1 b1,1 c2,2 a2,2 b2,2 c2,3 a2,3 b3,3 c2,4 a2,4 b4,4 + + + + c3,1 a3,2 b2,1 c3,2 a3,3 b3,2 c3,3 a3,4 b4,3 c3,4 a3,1 + + + 国家高性能计算中心（合肥） 2019/1/3

Systolic乘法 Step 7 b1,4 c1,1 c1,2 c1,3 a1,1 b1,3 c1,4 a1,2 b2,4 c2,1 + + c2,1 c2,2 a2,1 b1,2 c2,3 a2,2 b2,3 c2,4 a2,3 b3,4 + + + c3,1 a3,1 b1,1 c3,2 a3,2 b3,2 c3,3 a3,3 b3,3 c3,4 a3,4 b4,4 + + + + 国家高性能计算中心（合肥） 2019/1/3

Systolic乘法 Step 8 c1,1 c1,2 c1,3 c1,4 a1,1 b1,4 c2,1 c2,2 c2,3 a2,1 + c2,1 c2,2 c2,3 a2,1 b1,3 c2,4 a2,2 b2,4 + + c3,1 c3,2 a3,1 b1,2 c3,3 a3,2 b2,3 c3,4 a3,3 b3,4 + + + 国家高性能计算中心（合肥） 2019/1/3

Systolic乘法 Step 9 c1,1 c1,2 c1,3 c1,4 c2,1 c2,2 c2,3 c2,4 a2,1 b1,4 + c3,1 c3,2 c3,3 a3,1 b1,3 c3,4 a3,2 b2,4 + + 国家高性能计算中心（合肥） 2019/1/3

Systolic乘法 Step 10 c1,1 c1,2 c1,3 c1,4 c2,1 c2,2 c2,3 c2,4 c3,1 c3,2 a3,1 b1,4 + 国家高性能计算中心（合肥） 2019/1/3

Systolic乘法 Over c1,1 = a1,1 b1,1 + a1,2 b2,1 + a1,3 b3,1 + a1,4 b4,1 ………… c3,4 = a3,1 b1,4 + a3,2 b2,4 + a3,3 b3,4 + a3,4 b4,4 Over 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 国家高性能计算中心（合肥） 2019/1/3

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 End 国家高性能计算中心（合肥） 2019/1/3

9. 4 矩阵乘法 9. 4. 1 简单并行分块乘法 9. 4. 2 Cannon乘法 9. 4. 3 Fox乘法 9. 4 9.4 矩阵乘法 9.4.1 简单并行分块乘法 9.4.2 Cannon乘法 9.4.3 Fox乘法 9.4.4 Systolic乘法 9.4.5 DNS乘法

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方向进行单点积累求和; 国家高性能计算中心（合肥） 2019/1/3

示例 k j i C00=1×(-5)+2×7=9 C01=1×(-6)+2×8=10 C10=3×(-5)+4×7=13 A 101 111 P5 P7 100 110 P4 P6 001 011 P1 P3 000 010 P0 P2 A 初始加载 (b)A,B沿k维复制 (c)A沿j维复制 k j i B B (d)B沿i维复制 (e)点积 (f)沿k维求和 国家高性能计算中心（合肥） 2019/1/3

DNS乘法 算法描述: //令r(m)表示r的第m位取反； //{p, rm=d}表示r(0≤r≤p-1)的集合, 这里r的二 //输入: 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 (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) (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 (5)for m=2q to 3q-1 do //求和,m=2 for r=0 to p-1 par-do Cr=Cr+Cr(m) End 国家高性能计算中心（合肥） 2019/1/3

DNS乘法 示例 国家高性能计算中心（合肥） 2019/1/3

DNS乘法 k j i 国家高性能计算中心（合肥） 2019/1/3