并行计算 Parallel Computing 主讲人 孙广中 Spring, 2018

第三篇 并行数值算法 第九章 稠密矩阵运算 第十章 线性方程组的求解 第十一章 快速傅里叶变换 第十二章 数值计算的基本支撑技术 第三篇 并行数值算法 第九章 稠密矩阵运算 第十章 线性方程组的求解 第十一章 快速傅里叶变换 第十二章 数值计算的基本支撑技术

第九章 稠密矩阵运算 9. 1 矩阵的划分 9. 1. 1 带状划分 9. 1. 2 棋盘划分 9. 2 矩阵转置 9 第九章 稠密矩阵运算 9.1 矩阵的划分 9.1.1 带状划分 9.1.2 棋盘划分 9.2 矩阵转置 9.3 矩阵-向量乘法 9.4 矩阵乘法

划分方法 带状划分(striped partitioning): one dimensional, row or column, block or cyclic 棋盘划分(checkerboard partitioning): two dimensional, block or cyclic 国家高性能计算中心（合肥）

第九章 稠密矩阵运算 9. 1 矩阵的划分 9. 1. 1 带状划分 9. 1. 2 棋盘划分 9. 2 矩阵转置 9 第九章 稠密矩阵运算 9.1 矩阵的划分 9.1.1 带状划分 9.1.2 棋盘划分 9.2 矩阵转置 9.3 矩阵-向量乘法 9.4 矩阵乘法

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

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

第九章 稠密矩阵运算 9. 1 矩阵的划分 9. 1. 1 带状划分 9. 1. 2 棋盘划分 9. 2 矩阵转置 9 第九章 稠密矩阵运算 9.1 矩阵的划分 9.1.1 带状划分 9.1.2 棋盘划分 9.2 矩阵转置 9.3 矩阵-向量乘法 9.4 矩阵乘法

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

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

第九章 稠密矩阵运算 9. 1 矩阵的划分 9. 2 矩阵转置 9. 2. 1 棋盘划分的矩阵转置 9. 2. 2 带状划分的矩阵转置 9 第九章 稠密矩阵运算 9.1 矩阵的划分 9.2 矩阵转置 9.2.1 棋盘划分的矩阵转置 9.2.2 带状划分的矩阵转置 9.3 矩阵-向量乘法 9.4 矩阵乘法

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

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

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

棋盘划分的矩阵转置（4） 超立方连接：示例 3 6 7 4 5 2 1 14 15 12 13 10 11 8 9 P (0,5) (1,5) P (2,5) (3,5) (5,5) (4,5) (7,5) (6,5) (0,0) (0,1) (1,0) (1,1) (0,2) (0,3) (1,2) (1,3) (0,4) (1,4) (0,6) (0,7) (1,6) (1,7) (2,0) (2,1) (3,0) (3,1) (2,2) (2,3) (3,2) (3,3) (2,4) (3,4) (2,6) (2,7) (3,6) (3,7) (5,0) (5,1) (4,0) (4,1) (5,2) (5,3) (4,2) (4,3) (5,4) (4,4) (5,6) (5,7) (4,6) (4,7) (7,0) (7,1) (6,0) (6,1) (7,2) (7,3) (6,2) (6,3) (7,4) (6,4) (7,6) (7,7) (6,6) (6,7) (b) (a) 4 8 12 1 5 9 2 6 10 14 3 7 11 15 13 3 6 7 4 5 2 1 14 15 12 13 10 11 8 9 国家高性能计算中心（合肥）

第九章 稠密矩阵运算 9. 1 矩阵的划分 9. 2 矩阵转置 9. 2. 1 棋盘划分的矩阵转置 9. 2. 2 带状划分的矩阵转置 9 第九章 稠密矩阵运算 9.1 矩阵的划分 9.2 矩阵转置 9.2.1 棋盘划分的矩阵转置 9.2.2 带状划分的矩阵转置 9.3 矩阵-向量乘法 9.4 矩阵乘法

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

第九章 稠密矩阵运算 9. 1 矩阵的划分 9. 2 矩阵转置 9. 3 矩阵-向量乘法 9. 3. 1 带状划分的矩阵-向量乘法 9. 3 第九章 稠密矩阵运算 9.1 矩阵的划分 9.2 矩阵转置 9.3 矩阵-向量乘法 9.3.1 带状划分的矩阵-向量乘法 9.3.2 棋盘划分的矩阵-向量乘法 9.3.3 矩阵-向量的脉动乘法 9.4 矩阵乘法

矩阵-向量乘法 求Y=AX 串行算法计算时间t(n)=O(n2) 国家高性能计算中心（合肥）

第九章 稠密矩阵运算 9. 1 矩阵的划分 9. 2 矩阵转置 9. 3 矩阵-向量乘法 9. 3. 1 带状划分的矩阵-向量乘法 9. 3 第九章 稠密矩阵运算 9.1 矩阵的划分 9.2 矩阵转置 9.3 矩阵-向量乘法 9.3.1 带状划分的矩阵-向量乘法 9.3.2 棋盘划分的矩阵-向量乘法 9.3.3 矩阵-向量的脉动乘法 9.4 矩阵乘法

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

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

第九章 稠密矩阵运算 9. 1 矩阵的划分 9. 2 矩阵转置 9. 3 矩阵-向量乘法 9. 3. 1 带状划分的矩阵-向量乘法 9. 3 第九章 稠密矩阵运算 9.1 矩阵的划分 9.2 矩阵转置 9.3 矩阵-向量乘法 9.3.1 带状划分的矩阵-向量乘法 9.3.2 棋盘划分的矩阵-向量乘法 9.3.3 矩阵-向量的脉动乘法 9.4 矩阵乘法

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

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

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

第九章 稠密矩阵运算 9. 1 矩阵的划分 9. 2 矩阵转置 9. 3 矩阵-向量乘法 9. 3. 1 带状划分的矩阵-向量乘法 9. 3 第九章 稠密矩阵运算 9.1 矩阵的划分 9.2 矩阵转置 9.3 矩阵-向量乘法 9.3.1 带状划分的矩阵-向量乘法 9.3.2 棋盘划分的矩阵-向量乘法 9.3.3 矩阵-向量的脉动乘法 9.4 矩阵乘法

矩阵-向量乘法的脉动算法（1） 示例 国家高性能计算中心（合肥）

矩阵-向量乘法的脉动算法（2） 示例 国家高性能计算中心（合肥）

第九章 稠密矩阵运算 9. 1 矩阵的划分 9. 2 矩阵转置 9. 3 矩阵-向量乘法 9. 4 矩阵乘法 9. 4 第九章 稠密矩阵运算 9.1 矩阵的划分 9.2 矩阵转置 9.3 矩阵-向量乘法 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中元素的第2下标与B中元素的第1下标相一致（对准） 国家高性能计算中心（合肥）

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

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

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

第九章 稠密矩阵运算 9. 1 矩阵的划分 9. 2 矩阵转置 9. 3 矩阵-向量乘法 9. 4 矩阵乘法 9. 4 第九章 稠密矩阵运算 9.1 矩阵的划分 9.2 矩阵转置 9.3 矩阵-向量乘法 9.4 矩阵乘法 9.4.1 简单并行分块乘法 9.4.2 Cannon乘法 9.4.3 Fox乘法 9.4.4 Systolic乘法 9.4.5 DNS乘法

Cannon乘法（1） 分块: 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 国家高性能计算中心（合肥）

Cannon乘法（2） 算法原理 (非形式描述，1969年) 所有块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的每个块向上循环移动一步; ④转②执行 次; 国家高性能计算中心（合肥）

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

Cannon乘法（3） 示例： 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 国家高性能计算中心（合肥）

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

Cannon乘法（5） 算法描述: 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 初步的时间分析： 国家高性能计算中心（合肥）

Cannon乘法（6） 时间分析 (1)超立方连接: ②和③执行 次，所以运行时间为 (2)二维网孔连接，CT选路模式: ②和③执行 次，所以运行时间为 (2)二维网孔连接，CT选路模式: 国家高性能计算中心（合肥）

第九章 稠密矩阵运算 9. 1 矩阵的划分 9. 2 矩阵转置 9. 3 矩阵-向量乘法 9. 4 矩阵乘法 9. 4 第九章 稠密矩阵运算 9.1 矩阵的划分 9.2 矩阵转置 9.3 矩阵-向量乘法 9.4 矩阵乘法 9.4.1 简单并行分块乘法 9.4.2 Cannon乘法 9.4.3 Fox乘法 9.4.4 Systolic乘法 9.4.5 DNS乘法

Fox乘法（1） 分块: 同Cannon分块算法 算法原理（1987年） ①Ai,i向所在行的其他处理器 进行一到多播送； 有的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 国家高性能计算中心（合肥）

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

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

Fox乘法（4） 运行时间 (1)超立方连接: ②、③和④(含①)执行 次，所以运行时间为 当p=n2时, t(n)=O(nlogn) ②、③和④(含①)执行 次，所以运行时间为 当p=n2时, t(n)=O(nlogn) (2)二维网孔连接，CT选路模式(思考?) 国家高性能计算中心（合肥）

第九章 稠密矩阵运算 9. 1 矩阵的划分 9. 2 矩阵转置 9. 3 矩阵-向量乘法 9. 4 矩阵乘法 9. 4 第九章 稠密矩阵运算 9.1 矩阵的划分 9.2 矩阵转置 9.3 矩阵-向量乘法 9.4 矩阵乘法 9.4.1 简单并行分块乘法 9.4.2 Cannon乘法 9.4.3 Fox乘法 9.4.4 Systolic乘法 9.4.5 DNS乘法

Systolic乘法（1） 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 国家高性能计算中心（合肥）

Systolic乘法（2） 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 国家高性能计算中心（合肥）

Systolic乘法（3） 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 a1,4 c1,2 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 国家高性能计算中心（合肥）

Systolic乘法（4） 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 a2,4 c2,2 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 + 国家高性能计算中心（合肥）

Systolic乘法（5） Step 5 b1,4 b1,3 b2,4 b1,2 b2,3 b3,4 c1,1 a1,1 b1,1 c1,2 + + + + c2,1 a2,2 b2,1 a2,3 c2,2 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 + + 国家高性能计算中心（合肥）

Systolic乘法（6） Step 6 b1,4 b1,3 b2,4 c1,1 a1,1 c1,2 b1,2 c1,3 a1,2 b2,3 + + + c2,1 a2,1 b1,1 a2,2 c2,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 + + + 国家高性能计算中心（合肥）

Systolic乘法（7） Step 7 b1,4 c1,1 c1,2 c1,3 a1,1 b1,3 c1,4 a1,2 b2,4 c2,1 + + c2,1 a2,1 c2,2 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 + + + + 国家高性能计算中心（合肥）

Systolic乘法（8） 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 + + + 国家高性能计算中心（合肥）

Systolic乘法（9） 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 + + 国家高性能计算中心（合肥）

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

Systolic乘法（11） 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 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 国家高性能计算中心（合肥）

Systolic乘法（12） Systolic算法（H.T. Kung） //输入: 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 国家高性能计算中心（合肥）

第九章 稠密矩阵运算 9. 1 矩阵的划分 9. 2 矩阵转置 9. 3 矩阵-向量乘法 9. 4 矩阵乘法 9. 4 第九章 稠密矩阵运算 9.1 矩阵的划分 9.2 矩阵转置 9.3 矩阵-向量乘法 9.4 矩阵乘法 9.4.1 简单并行分块乘法 9.4.2 Cannon乘法 9.4.3 Fox乘法 9.4.4 Systolic乘法 9.4.5 DNS乘法

DNS乘法（1） Motivation: From a good and common idea j C A B i +  ... = 1 2 3 n ... = 国家高性能计算中心（合肥）

DNS乘法（2） Motivation: From a good and common idea(Cont.) How to use processors more effectively and practically? x 1 x 2 x 3 x n ... + + n processors for each result element implies n x n2 = n3 ... = + time is log2 n 国家高性能计算中心（合肥）

DNS乘法（3） 背景: 1981年由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方向进行单点积累求和; 国家高性能计算中心（合肥）

示例 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维求和 国家高性能计算中心（合肥）

DNS乘法（5） 算法描述: //令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 国家高性能计算中心（合肥）

DNS乘法（6） 示例 国家高性能计算中心（合肥）

DNS乘法（7） k j i 国家高性能计算中心（合肥）