并行计算 课件.ppt

并行计算Parallel Computing主讲人 孙广中Spring,2018Spring,2018第三篇 并行数值算法 第九章 稠密矩阵运算 第十章 线性方程组的求解 第十一章 快速傅里叶变换 第十二章 数值计算的基本支撑技术 第九章 稠密矩阵运算 9.1 矩阵的划分 9.1.1 带状划分 9.1.2 棋盘划分 9.2 矩阵转置 9.3 矩阵-向量乘法 9.4 矩阵乘法国家高性能计算中心（合肥）4 划分方法带状划分带状划分(striped partitioning):(striped partitioning):one dimensional,row or column,one dimensional,row or column,block or cyclic block or cyclic棋盘划分棋盘划分(checkerboard partitioning):(checkerboard partitioning):two dimensional,block or cyclic two dimensional,block or cyclic第九章 稠密矩阵运算 9.1 矩阵的划分 9.1.1 带状划分 9.1.2 棋盘划分 9.2 矩阵转置 9.3 矩阵-向量乘法 9.4 矩阵乘法国家高性能计算中心（合肥）6 带状划分（1）16161616阶矩阵，阶矩阵，p=4 p=4 列块带状划分列块带状划分 行循环带状划分行循环带状划分国家高性能计算中心（合肥）7 带状划分（2）示例：示例：p p3 3，27 2727 27矩阵的矩阵的3 3种带状划分种带状划分 第九章 稠密矩阵运算 9.1 矩阵的划分 9.1.1 带状划分 9.1.2 棋盘划分 9.2 矩阵转置 9.3 矩阵-向量乘法 9.4 矩阵乘法国家高性能计算中心（合肥）9 棋盘划分（1）8888阶矩阵，阶矩阵，p=16 p=16 块棋盘划分块棋盘划分 循环棋盘划分循环棋盘划分国家高性能计算中心（合肥）10 棋盘划分（2）示例：示例：p p4 4，16161616矩阵的矩阵的3 3种棋盘划分种棋盘划分 第九章 稠密矩阵运算 9.1 矩阵的划分 9.2 矩阵转置 9.2.1 棋盘划分的矩阵转置 9.2.2 带状划分的矩阵转置 9.3 矩阵-向量乘法 9.4 矩阵乘法国家高性能计算中心（合肥）12 棋盘划分的矩阵转置（1）网孔连接 情形情形1:p=n1:p=n2 2。通讯步通讯步 转置后转置后国家高性能计算中心（合肥）13 棋盘划分的矩阵转置（2）情形情形2:pn2:pn2 2。-划分划分:-算法算法:按按meshmesh连接进行块转置连接进行块转置(不同处理器间不同处理器间)进行块内转置进行块内转置(同一处理器内同一处理器内)通讯步通讯步 转置后转置后国家高性能计算中心（合肥）14 棋盘划分的矩阵转置（3）超立方连接 划分划分:算法算法:对对A Aij ij递归应用递归应用进行转置，直至分块矩阵的元素处于同一处理器；进行转置，直至分块矩阵的元素处于同一处理器；进行同一处理器的内部转置。进行同一处理器的内部转置。运行时间运行时间:国家高性能计算中心（合肥）15 棋盘划分的矩阵转置（4）超立方连接：示例示例 3674520114151213101189(0,5)(1,5)P PP PP P(2,5)(3,5)P PP PP PP P(5,5)(4,5)P PP PP PP P(7,5)(6,5)P P(0,0)(0,1)(1,0)(1,1)(0,2)(0,3)(1,2)(1,3)P P(0,4)(0,5)(1,4)(1,5)P P(0,6)(0,7)(1,6)(1,7)P P(2,0)(2,1)(3,0)(3,1)P P(2,2)(2,3)(3,2)(3,3)P P(2,4)(2,5)(3,4)(3,5)P 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 P(5,4)(5,5)(4,4)(4,5)P P(5,6)(5,7)(4,6)(4,7)P P(7,0)(7,1)(6,0)(6,1)P P(7,2)(7,3)(6,2)(6,3)P P(7,4)(7,5)(6,4)(6,5)P P(7,6)(7,7)(6,6)(6,7)(b)(a)048121592610143711150841219513621014731115(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 PP PP PP PP PP PP PP 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)13(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.1 矩阵的划分 9.2 矩阵转置 9.2.1 棋盘划分的矩阵转置 9.2.2 带状划分的矩阵转置 9.3 矩阵-向量乘法 9.4 矩阵乘法国家高性能计算中心（合肥）17 带状划分的矩阵转置 划分划分:A:Annnn分成分成p p个个(n/p)n(n/p)n大小的带大小的带 算法算法:PiPi有有p-1p-1个个(n/p)(n/p)(n/p)(n/p)大小子块发送到另外大小子块发送到另外p-1p-1个处理器中个处理器中;每个处理器本地交换相应的元素；每个处理器本地交换相应的元素；时间分析？时间分析？第九章 稠密矩阵运算 9.1 矩阵的划分 9.2 矩阵转置 9.3 矩阵-向量乘法 9.3.1 带状划分的矩阵-向量乘法 9.3.2 棋盘划分的矩阵-向量乘法 9.3.3 矩阵-向量的脉动乘法 9.4 矩阵乘法国家高性能计算中心（合肥）19矩阵-向量乘法 求求Y=AXY=AX 串行算法计算时间串行算法计算时间t(n)=O(nt(n)=O(n2 2)第九章 稠密矩阵运算 9.1 矩阵的划分 9.2 矩阵转置 9.3 矩阵-向量乘法 9.3.1 带状划分的矩阵-向量乘法 9.3.2 棋盘划分的矩阵-向量乘法 9.3.3 矩阵-向量的脉动乘法 9.4 矩阵乘法国家高性能计算中心（合肥）21 带状划分的矩阵-向量乘法（1）划分划分(行带状划分行带状划分):P):Pi i存放存放x xi i和和a ai,0i,0,a,ai,1i,1,a,ai,n-1i,n-1,并输出并输出y yi i 算法算法:对对p=np=n情形情形 每个每个P Pi i向其他处理器播送向其他处理器播送x xi i(多到多播送多到多播送)；每个每个P Pi i做相应计算；做相应计算；注注:对对pnpn情形情形,算法中算法中P Pi i要播送要播送X X中相应的中相应的n/pn/p个分量个分量 (1)(1)超立方连接的计算时间超立方连接的计算时间 (2)(2)网孔连接的计算时间网孔连接的计算时间国家高性能计算中心（合肥）22 带状划分的矩阵-向量乘法（2）示例示例第九章 稠密矩阵运算 9.1 矩阵的划分 9.2 矩阵转置 9.3 矩阵-向量乘法 9.3.1 带状划分的矩阵-向量乘法 9.3.2 棋盘划分的矩阵-向量乘法 9.3.3 矩阵-向量的脉动乘法 9.4 矩阵乘法国家高性能计算中心（合肥）24 棋盘划分的矩阵-向量乘法（1）划分划分(块棋盘划分块棋盘划分):):P Pij ij存放存放a ai,ji,j,x,xi i置入置入P Pi,ii,i中中算法算法:对对p=np=n2 2情形情形 每个每个P Pi,ii,i向向P Pj,ij,i播送播送x xi i(一到多播送一到多播送)；按行方向进行乘按行方向进行乘-加与积累运算，最后一列加与积累运算，最后一列P Pi,n-1i,n-1收集的结收集的结果为果为y yi i；注注:对对pnp p)(n p)P0,0P1,0P2,0P3,0P0,1P1,1P2,1P3,1P0,2P1,2P2,2P3,2P0,3P1,3P2,3P3,3国家高性能计算中心（合肥）37Cannon乘法（2）算法原理(非形式描述，非形式描述，19691969年年)所有块所有块A Ai,ji,j(0i,j )(0i,j )向左循环移动向左循环移动i i步步(按行移位按行移位)；所有块所有块B Bi,ji,j(0i,j )(0i,j )向上循环移动向上循环移动j j步步(按列移位按列移位)；所有处理器所有处理器P Pi,ji,j做执行做执行A Ai,ji,j和和B Bi,ji,j的乘的乘-加运算加运算;AA的每个块向左循环移动一步的每个块向左循环移动一步;B B的每个块向上循环移动一步的每个块向上循环移动一步;转转执行执行 次次;国家高性能计算中心（合肥）38Cannon乘法（2）示例：A A4444，B B4444,p=16,p=16 A0,0A1,0A2,0A3,0A0,1A1,1A2,1A3,1A0,2A1,2A2,2A3,2A0,3A1,3A2,3A3,3B0,0B1,0B2,0B3,0B0,1B1,1B2,1B3,1B0,2B1,2B2,2B3,2B0,3B1,3B2,3B3,3Initial alignment of AInitial alignment of B国家高性能计算中心（合肥）39Cannon乘法（3）示例：A A4444，B B4444,p=16,p=16 A and B after initial alignment and shifts after every stepA0,0B0,0A1,1B1,0A2,2B2,0A3,3B3,0A0,1B1,1A1,2B2,1A2,3B3,1A3,0B0,1A0,2B2,2A1,3B3,2A2,0B0,2A3,1B1,2A0,3B3,3A1,0B0,3A2,1B1,3A3,2B2,3国家高性能计算中心（合肥）40Cannon乘法（4）示例：A A4444，B B4444,p=16,p=16 After first shiftA0,1B1,0A1,2B2,0A2,3B3,0A3,0B0,0A0,2B2,1A1,3B3,1A2,0B0,1A3,1B3,1A0,3B3,2A1,0B0,2A2,1B1,2A3,2B2,2A0,0B0,3A1,1B1,3A2,2B2,3A3,3B3,3After second shiftA0,2B2,0A1,3B3,0A2,0B0,0A3,1B1,0A0,3B3,1A1,0B0,1A2,1B1,1A3,2B2,1A0,0B0,2A1,1B1,2A2,2B2,2A3,3B3,2A0,1B1,3A1,2B2,3A2,3B3,3A3,0B0,3After third shiftA0,3B3,0A1,0B0,0A2,1B1,0A3,2B2,0A0,0B0,1A1,1B1,1A2,2B2,1A3,3B3,1A0,1B1,2A1,2B2,2A2,3B3,2A3,0B0,2A0,2B2,3A1,3B3,3A2,0B0,3A3,1B1,3国家高性能计算中心（合肥）41Cannon乘法（5）算法描述:Cannon分块乘法算法 /输入输入:A:An nn n,B,Bn nn n;输出输出:C:Cn nn n Begin Begin (1)for k=0 to do (1)for k=0 to do for all P for all Pi,ji,j par-do par-do (i)if ik then (i)if ik then A Ai,ji,j A Ai,(j+1)modi,(j+1)mod endif endif (ii)if jk then (ii)if jk then B Bi,ji,j B B(i+1)mod ,j(i+1)mod ,j endif endif endfor endfor endfor endfor (2)for all P (2)for all Pi,ji,j par-do C par-do Ci,ji,j=0 endfor =0 endfor (3)for k=0 to do(3)for k=0 to do for all P for all Pi,ji,j par-do par-do (i)C (i)Ci,ji,j=C=Ci,ji,j+A+Ai,ji,jB Bi,ji,j (ii)A (ii)Ai,ji,j A Ai,(j+1)modi,(j+1)mod (iii)B (iii)Bi,ji,j B B(i+1)mod ,j(i+1)mod ,j endfor endfor endfor endfor End End初步的时间分析：初步的时间分析：国家高性能计算中心（合肥）42Cannon乘法（6）时间分析 (1)(1)超立方连接超立方连接:和和执行执行 次，所以运行时间为次，所以运行时间为 (2)(2)二维网孔连接，二维网孔连接，CTCT选路模式选路模式:和和执行执行 次，所以运行时间为次，所以运行时间为 第九章 稠密矩阵运算 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乘法国家高性能计算中心（合肥）44Fox乘法（1）分块:同同CannonCannon分块算法分块算法算法原理（19871987年）年）A Ai,ii,i向所在行的其他处理器向所在行的其他处理器进行一到多播送；进行一到多播送；各处理器将收到的各处理器将收到的A A块与原块与原有的有的B B块进行乘块进行乘-加运算加运算;BB块向上循环移动一步块向上循环移动一步;如果如果A Ai,ji,j是上次第是上次第i i行播送的块，本次选择行播送的块，本次选择 向所向所 在行的其他处理器进行一到多播送在行的其他处理器进行一到多播送;转转执行执行 次次;A0,0B0,0A1,0B1,0A2,0B2,0A3,0B3,0A0,1B0,1A1,1B1,1A2,1B2,1A3,1B3,1A0,2B0,2A1,2B1,2A2,2B2,2A3,2B3,2A0,3B0,3A1,3B1,3A2,3B2,3A3,3B3,3国家高性能计算中心（合肥）45Fox乘法（2）示例：A A4444，B B4444,p=16,p=16 (a)(b)(a)(b)A0,0B0,0B1,0B2,0B3,0B0,1A1,1B1,1B2,1B3,1B0,2B1,2A2,2B2,2B3,2B0,3B1,3B2,3A3,3B3,3B1,0B2,0B3,0A0,1B1,1B2,1B3,1B0,1B1,2B3,2B0,2B1,3B2,3B0,3A1,2B2,2A2,3B3,3A3,0B0,0国家高性能计算中心（合肥）46Fox乘法（3）示例：A A4444，B B4444,p=16,p=16 (c)(d)(c)(d)B2,0B3,0B2,1B3,1B0,1B3,2B0,2B1,2B2,3B0,3B1,3B3,0B1,0B3,1B0,1B2,1B3,2B1,2B0,3B2,3B0,2B1,3B2,0A0,2B2,2A1,3B3,3A2,0B0,0B1,0A3,1B1,1A0,3B3,3A1,0B0,2A2,1B1,1A3,2B2,2国家高性能计算中心（合肥）47Fox乘法（4）运行时间运行时间 (1)(1)超立方连接超立方连接:、和和(含含)执行执行 次，所以运行时间为次，所以运行时间为 当当p=np=n2 2时时,t(n)=O(nlogn),t(n)=O(nlogn)(2)(2)二维网孔连接，二维网孔连接，CTCT选路模式选路模式(思考思考?)?)第九章 稠密矩阵运算 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乘法国家高性能计算中心（合肥）49Systolic乘法（1）a1,4b4,1b3,1b2,1b2,2b4,2b3,2b2,3b3,3b4,3b2,4b3,4b4,4a1,3a1,1a1,2a2,4a2,1a2,2a2,3a3,1a3,2a3,3a3,4b1,1b1,2b1,3b1,4Step 1P1,1c1,1P1,2c1,2P1,3c1,3P1,4c1,4P2,1c2,1P2,2c2,2P2,3c2,3P2,4c2,4P3,1c3,1P3,2c3,2P3,3c3,3P3,4c3,4国家高性能计算中心（合肥）50Systolic乘法（2）c1,1c1,2c1,3c1,4c2,1c2,2c2,3c2,4c3,1c3,2c3,3c3,4b3,1b2,1b2,2b4,2b3,2b2,3b3,3b4,3b2,4b3,4b4,4a1,3a1,1a1,2a2,4a2,1a2,2a2,3a3,1a3,2a3,3a3,4b1,1b1,2b1,3b1,4a1,4b4,1+Step 2国家高性能计算中心（合肥）51Systolic乘法（3）c1,1c1,2c1,3c1,4c2,1c2,2c2,3c2,4c3,1c3,2c3,3c3,4b2,1b2,2b3,2b2,3b3,3b4,3b2,4b3,4b4,4a1,1a1,2a2,1a2,2a2,3a3,1a3,2a3,3a3,4b1,1b1,2b1,3b1,4a1,3b3,1+a1,4b4,2+a2,4b4,1+Step 3国家高性能计算中心（合肥）52Systolic乘法（4）c1,1c1,2c1,3c1,4c2,1c2,2c2,3c2,4c3,1c3,2c3,3c3,4b2,2b2,3b3,3b2,4b3,4b4,4a1,1a2,1a2,2a3,1a3,2a3,3b1,1b1,2b1,3b1,4a1,2b2,1+a1,3b3,2+a2,3b3,1+a1,4b4,3+a3,4b4,1+a2,4b4,2+Step 4国家高性能计算中心（合肥）53Systolic乘法（5）c1,1c1,2c1,3c1,4c2,1c2,2c2,3c2,4c3,1c3,2c3,3c3,4b2,3b2,4b3,4a2,1a3,1a3,2b1,2b1,3b1,4a1,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国家高性能计算中心（合肥）54Systolic乘法（6）c1,1c1,2c1,3c1,4c2,1c2,2c2,3c2,4c3,1c3,2c3,3c3,4b2,4a3,1b1,3b1,4a1,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国家高性能计算中心（合肥）55Systolic乘法（7）c1,1c1,2c1,3c1,4c2,1c2,2c2,3c2,4c3,1c3,2c3,3c3,4b1,4a1,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国家高性能计算中心（合肥）56Systolic乘法（8）c1,1c1,2c1,3c1,4c2,1c2,2c2,3c2,4c3,1c3,2c3,3c3,4a1,1b1,4+a2,1b1,3+a3,1b1,2+a2,2b2,4+a3,2b2,3+a3,3b3,4+Step 8国家高性能计算中心（合肥）57Systolic乘法（9）c1,1c1,2c1,3c1,4c2,1c2,2c2,3c2,4c3,1c3,2c3,3c3,4a2,1b1,4+a3,1b1,3+a3,2b2,4+Step 9国家高性能计算中心（合肥）58Systolic乘法（10）c1,1c1,2c1,3c1,4c2,1c2,2c2,3c2,4c3,1c3,2c3,3c3,4a3,1b1,4+Step 10国家高性能计算中心（合肥）59Systolic乘法（11）P1,1c1,1P1,2c1,2P1,3c1,3P1,4c1,4P2,1c2,1P2,2c2,2P2,3c2,3P2,4c2,4P3,1c3,1P3,2c3,2P3,3c3,3P3,4c3,4Overc1,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 国家高性能计算中心（合肥）60Systolic乘法（12）Systolic算法（H.T.KungH.T.Kung）/输入输入:A:Ammn n,B,Bn nk k;输出输出:C:Cmmk k Begin Begin for i=1 to m par-do for i=1 to m par-do for j=1 to k par-do for j=1 to k par-do (i)c (i)ci,ji,j=0=0 (ii)while P (ii)while Pi,j i,j 收到收到a a和和b b时时 dodo c ci,ji,j=c=ci,ji,j+ab+ab if i m then if i m then 发送发送b b给给P Pi+1,ji+1,j endif endif if j k then if j k then 发送发送a a给给P Pi,j+1i,j+1 endif endif endwhile endwhile endfor endfor endfor endfor End End第九章 稠密矩阵运算 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乘法国家高性能计算中心（合肥）62DNS乘法（1）Motivation:From a good and common ideaMotivation:From a good and common idea ji+11 22 33 nn+.+=ABC国家高性能计算中心（合肥）63DNS乘法（2）Motivation:From a good and common idea(Cont.)Motivation:From a good and common idea(Cont.)How to use processors more effectively and practically?How to use processors more effectively and practically?x11x22x33xnn .=.+n processors for each result element implies n x n2=n3time is log2 n国家高性能计算中心（合肥）64DNS乘法（3）背景背景:1981 1981年由年由DekelDekel、NassimiNassimi和和SahniSahni提出的提出的SIMD-CCSIMD-CC上的矩阵乘法上的矩阵乘法,处理处理 器数目为器数目为n n3 3,运行时间为运行时间为O(logn),O(logn),是一种速度很快的算法。是一种速度很快的算法。基本思想基本思想:通过一到一和一到多的播送办法，使得处理器通过一到一和一到多的播送办法，使得处理器(k,i,j)(k,i,j)拥有拥有a ai,ki,k和和b bk,jk,j,进行本地相乘进行本地相乘,再沿再沿k k方向进行单点积累求和方向进行单点积累求和,结果存储在处理器结果存储在处理器(0,i,j)(0,i,j)中。中。处理器编号处理器编号:处理器数处理器数p=np=n3 3=(2=(2q q)3 3=2=23 3q q,处理器处理器P Pr r位于位置位于位置(k,i,j),(k,i,j),这里这里r=knr=kn2 2+in+j,(0i,j,kn-1)+in+j,(0i,j,kn-1)。位于。位于(k,i,j)(k,i,j)的处理器的处理器P Pr r的三个寄存器的三个寄存器 A Ar r,B,Br r,C,Cr r分别表示为分别表示为Ak,i,j,Bk,i,jAk,i,j,Bk,i,j和和Ck,i,j,Ck,i,j,初始时均为初始时均为0 0。算法算法:初始时初始时a ai,ji,j和和b bi,ji,j存储于寄存器存储于寄存器A0,i,jA0,i,j和和B0,i,j;B0,i,j;数据复制数据复制:A,B:A,B同时在同时在k k维复制维复制(一到一播送一到一播送)；A A在在j j维复制维复制(一到多播送一到多播送);B);B在在i i维复制维复制(一到多播送一到多播送););相乘运算相乘运算:所有处理器的所有处理器的A A、B B寄存器两两相乘寄存器两两相乘;求和运算求和运算:沿沿k k方向进行单点积累求和方向进行单点积累求和;国家高性能计算中心（合肥）65示例示例 C00=1(-5)+27=9C00=1(-5)+27=9C01=1(-6)+28=10C01=1(-6)+28=10C10=3(-5)+47=13C10=3(-5)+47=13C11=3(-6)+48=14C11=3(-6)+48=14 kji初始加载(b)A,B沿k维复制(c)A沿j维复制(d)B沿i维复制(e)点积(f)沿k维求和BBAA000010001011100110101111P0P1P2P3P4P5P6P7国家高性能计算中心（合肥）66算法描述算法描述:/令令r r(m)(m)表示表示r r的第的第mm位取反；位取反；/p,rp,rmm=d=d表示表示r(0r(0rrp-1p-1)的集合的集合,这里这里r r的二的二 /进制第进制第mm位为位为d;d;/输入输入:A:Annnn,B,Bnnnn;输出输出:C:Cnnnn Begin /Begin /以以n=2,p=8=2n=2,p=8=23 3举例举例,q=1,r=(r,q=1,r=(r2 2r r1 1r r0 0)2 2 (1)for m=3q-1 to 2q do /(1)for m=3q-1 to 2q do /按按k k维复制维复制A,B,A,B,m=2m=2 for all r in p,r for all r in p,rmm=0 par-do=0 par-do/r r2 2=0=0的的r r (1.1)A (1.1)Ar r(m)(m)A Ar r /A(100)A(100)A(000)A(000)等等 (1.2)B(1.2)Br r(m)(m)B Br r /B(100)B(100)B(000)B(000)等等 endforendfor endfor endfor (2)for m=q-1 to 0 do /(2)for m=q-1 to 0 do /按按j j维复制维复制A,A,m=0m=0 for all r in p,r for all r in p,rmm=r=r2q+m2q+m par-do/par-do/r r0 0=r=r2 2的的r r A Ar r(m)(m)A Ar r /A(001)A(001)A(000),A(100)A(000),A(100)A(101)A(101)endfor endfor /A(011)A(011)A(010),A(110)A(010),A(110)A(111)A(111)endfor endfor(3)for m=2q-1 to q do /(3)for m=2q-1 to q do /按按i i维复制维复制B,B,m=1m=1 for all r in p,r for all r in p,rmm=r=rq+mq+m par-do par-do/r r1 1=r=r2 2的的r r B Br r(m)(m)B Br r /B(010)B(010)B(000),B(100)B(000),B(100)B(110)B(110)endfor endfor /B(011)B(011)B(001),B(101)B(001),B(101)B(111)B(111)endfor endfor (4)for r=0 to p-1 par-do /(4)for r=0 to p-1 par-do /相乘相乘,all Pall Pr r C Cr r=A=Ar rBBr r endforendfor (5)for m=2q to 3q-1 do /(5)for m=2q to 3q-1 do /求和求和,m=2m=2 for r=0 to p-1 par-dofor r=0 to p-1 par-do C Cr r=C=Cr r+C+Cr r(m)(m)endfor endfor endfor endfor End EndDNS乘法（5）国家高性能计算中心（合肥）67示例示例 DNS乘法（6）国家高性能计算中心（合肥）68kjiDNS乘法（7）