Cannon乘法的MPI实现解析

1、并行算法实践要求学生在KD60实验平台上设计并行算法并实现。实验平台由一块处理板、一块监控板和一块背板等组成。处理板逻辑结构如图1所示。处理板承载4个处理单元，每个处理单元包括一个龙芯3号四核CPU、2GB DDR2内存、RTL8110千兆以太网卡芯片、BIOS Flash、串口收发芯片以及电源变换电路等。四个龙芯3号处理器通过HT总线实现互连。监控电路检测4个处理单元的状态，并实现对其控制。 图1 处理板逻辑结构实验平台的系统软件以开源软件为主（见图2），具有兼容性强、易维护、易升级、易使用等特点。处理单元操作系统为Debian GNU/Linux无盘系统，采用稳定高效的内核。图2 软件系统

2、结构要求选修并行算法实践同学在下面表1中选一个题目，（1）阐述基本原理（包括对算法改进和优化方法）；(2)根据KD60实验平台设计实验方法和步骤（包括主要调试过程要求拷屏）。(3) 数据及结果分析：根据不同的实验内容，记录具体的实验数据或程序运行结果（要求拷屏）。实验数据量较大时，最好制成表格形式。附上程序功能、模块说明、完整源代码，源代码中尽量多带注释；(4)分析和总结：对程序与算法性能改进结论，总结和体会。表1 并行算法实践题目 序号题目名称基本方法和内容要求1LU分解的Open MP实现编写LU分解的Open MP程序2KMP算法的Open MP实现编写KMP算法的OpenMP程序3高斯

3、消元法解线性方程组的Open MP实现编写高斯消元法解线性方程组的OpenMP程序4高斯消元法解线性方程组的MPI实现编写高斯消元法解线性方程组的MPI程序5高斯-塞德尔迭代解线性方程组的MPI实现编写高斯-塞德尔迭代解线性方程组的MPI程序6Cannon乘法的MPI实现编写Cannon乘法的MPI程序7LU分解的MPI实现编写LU分解的MPI程序8随机串匹配算法的MPI实现编写随机串匹配算法的MPI程序9单源最短路径Dijkstra 算法的MPI实现编写单源最短路径Dijkstra算法的MPI程序10快速排序算法的MPI实现编写快速排序算法的MPI程序11KMP串匹配的MPI实现编写KMP串

4、匹配算法的MPI程序图2 软件系统结构Cannon乘法的MPI实现及性能分析摘要：cannon算法是矩阵的并行乘法，属于数值并行算法MPI编程实现一篇，其中关于数值并行算法MPI编程由于要处理的数据量巨大，程序循环次数多，对于串行而言，处理时间将非常长，将其并行化非常必要。本文将矩阵数据进行棋盘划分成多个子矩阵，再分别指派给多个处理器，使个处理器并行运算。关键字：cannon乘法 并行计算 数据划分一、Cannon乘法的MPI实现基本原理Cannon乘法属于数值并行算法MPI编程实现一篇，其中关于数值并行算法MPI编程由于要处理的数据量巨大，程序循环次数多，对于串行而言，处理时间将非常长，使其

5、并行化的一般方法有：1）数据相关分析 2）数据划分和处理器指派 3）循环重构对原有程序并行化，首先要分析计算程序中所有语句间的依赖关系，这称之为相关分析。本项目Cannon乘法的mpi实现，是矩阵运算，阶往往都很高，而且行列之间数据依赖关系也不强，所以就对矩阵进行划分，然后指派给不同的处理器进行处理。最常用的矩阵划分有带状划分和块状划分。1 带状划分方法带状划分又叫行列划分，就是将矩阵整行或整列地分成若干组，各组指派给一个处理器。也可以将若干行或列指派给一个处理器，而且这些行和列可以是连续的，也可以是等间距的，前者称为块带状的，后者称为循环带状的。2. 块状划分方法块状划分又叫棋盘划分，就是将

6、矩阵划分成若干个子矩阵，每个子矩阵指派给一个处理器，此时任意处理器均不含整行或整列。和带状划分类似，棋盘划分也可分为块棋盘划分和循环棋盘划分。棋盘划分比带状划分可开发更高的并行度，Cannon乘法的mpi实现也正是基于棋盘划分的并行实现。循环重构是指在数据分解之后，相应地将串行程序循环部分进行重构，以实现这种划分所确定的并行计算，主要方法有1）循环交换 2）拉伸法 3）分裂法 4）轮转法 5）并列法 在三种程序并行化的方法中，数据相关分析和循环重构目的都是挖掘语句间的并行性，而数据划分和处理器指派则重在策略，宏观上挖掘并行性。Cannon算法是一种存储有效的算法，设矩阵和相乘。为了使两矩阵下标

7、满足相乘的要求，和带状的并行分块乘法不同，不是仅仅让B矩阵的各列块循环移动，而是有目的地让A的各行块以及B的各列块皆施行循环移位，从而实现对C的子块的计算。将矩阵A和B分成p个方块Aij和Bij，每块大小为，并将它们分配给个处理器。开始时处理器Pij存放块Aij和Bij，并负责计算块Cij，然后算法开始执行： 将块Aij向左循环移动i步；将块Bij向上循环移动j步；Pij执行乘加运算后将块Aij向左循环移动1步，块Bij向上循环移动1步；重复第步，总共执行次乘加运算和次块Aij和Bij的循环单步移位。二、Cannon乘法的MPI实现内容和步骤实验涉及内容主要有：1）数据划分和指派处理器最常用的

8、矩阵数据划分有带状划分和块状划分。棋盘划分比带状划分可开发更高的并行度，Cannon乘法的mpi实现也正是基于棋盘划分的并行实现。设有P个处理器，将矩阵A和B分成p个方块Aij和Bij，每块大小为，并将它们分配给个处理器。2） 子矩阵的循环移动处理器Pij存放块Aij和Bij，并负责计算块Cij，在使A矩阵的左右循环移动和B矩阵的上下循环移动时，为了避免在通信过程中发生死锁，奇数号及偶数号处理器的收发顺序被错开，使偶数号处理器先发送后接收；而奇数号处理器先将子矩阵块存于缓冲区Buffer中，然后接收编号在其后面的处理器所发送的子矩阵块，最后再将缓冲区中子矩阵块发送给编号在其前面的处理器。基本算

9、法如下：Begin(1)if (j=0) then /*最左端的子块*/(1.1)将所存的A的子块发送到同行最右端子块所在的处理器中(1.2)接收其右邻处理器中发来的A的子块end if (2)if (j = sqrt(p)-1) and (j mod 2 = 0) then /*最右端子块处理器且块列号为偶数*/(2.1)将所存的A的子块发送到其左邻处理器中(2.2)接收其同行最左端子块所在的处理器发来的A的子块end if (3)if (j = sqrt(p)-1) and (j mod 2 0) then /*最右端子块处理器且块列号为奇数*/(3.1)将所存的A的子块在缓冲区buffe

10、r中做备份(3.2)接收其同行最左端子块所在的处理器发来的A的子块(3.3)将在缓冲区buffer中所存的A的子块发送到其左邻处理器中end if (4)if (j sqrt(p)-1) and (j mod 2 = 0) and (j 0) then /*其余的偶数号处理器*/(4.1)将所存的A的子块发送到其左邻处理器中(4.2)接收其右邻处理器中发来的A的子块end if (5)if (j sqrt(p)-1) and (j mod 2 = 1) and (j 0) then /*其余的奇数号处理器*/(5.1)将所存的A的子块在缓冲区buffer中做备份(5.2)接收其右邻处理器中发来

11、的A的子块(5.3)将在缓冲区buffer中所存的A的子块发送到其左邻处理器中end if End实验步骤1） 登陆KD-60图 2.1 KD-60登陆界面2） 转至node80节点，上传程序输入命令：ssh loongsonnode80 和密码 进入图界面图2.2 转到节点80的界面再命令vim，进入vim编辑器加入程序,保存为cannon.c 3） 编译程序输入命令：mpicc cannon.c o cannon lm 在目录中查看，已成功。如下图图2.3 将程序保存并编译后界面4） 运行程序输入：mpirun np 4 cannon 4 ,其中第一个4是指定的处理器个数，第二个4是产生随

12、机矩阵的维数，这两个参数在实验过程中可以调整，但要求第一个参数即处理器的个数必须是一个数的平方数。输出：图2.4 cannon乘法运行结果图 2.4 并行程序运行界面 两个参数都是4， 分别输出两个随机矩阵和矩阵的乘积三、数据及结果1.下面列出了两组数据，分别是用一个处理器进行串行运算和四个处理器进行并行运算矩阵维数为200的计算时间比较。四个处理器处理阶数为200的矩阵相乘时，所花时间为：1.705844秒。单个处理器处理阶数为200的矩阵相乘时，所花时间为：3.727210秒。如图3.1 和图3.2所示。图3.1 四个处理器并行执行结果图图3.2 单个处理器串行执行结果图附：1. 程序模块

13、伪代码：输入：An×n，Bn×n输出：Cn×nBegin对所有处理器my_rank(my_rank=0,p-1)同时执行如下的算法:(1)计算子块的行号 i=my_rank/sqrt(p)计算子块的列号 j=my_rank mod sqrt(p)(2)for k=0 to -1 doif (i>k) then Leftmoveonestep(a) end if /* a循环左移至同行相邻处理器中*/if (j>k) then Upmoveonestep(b) end if /* b循环上移至同列相邻处理器中*/ end for (3)for i=0 t

14、o m-1 dofor j=0 to m-1 doci,j=0end forend for (4)for k=0 to -1 dofor i=0 to m-1 dofor j=0 to m-1 dofor k1=0 to m-1 doci,j= ci,j+ ai,k1* bk1,jend forend forend for Leftmoveonestep(a) /*子块a循环左移至同行相邻的处理器中*/Upmoveonestep(b) /*子块b循环上移至同列相邻的处理器中*/end for EndLeftmoveonestep(a) 见实验内容处附2.程序源码#include <std

15、lib.h>#include <string.h>#include <mpi.h>#include <time.h>#include <stdio.h>#include <math.h>/* 全局变量声明 */float *A, *B, *C; /* 总矩阵,C = A * B */float *a, *b, *c, *tmp_a, *tmp_b; /* a、b、c表分块，tmp_a、tmp_b表缓冲区 */int dg, dl, dl2,p, sp; /* dg:总矩阵维数;dl:矩阵块维数;dl2=dl*dl;p:处理器个数

16、;spsqrt(p) */int my_rank, my_row, my_col; /* my_rank:处理器ID;(my_row,my_col):处理器逻辑阵列坐标 */MPI_Status status;float starttime;float time1;/* *函数名: get_index *功能：处理器逻辑阵列坐标至rank号的转换 *输入：坐标、逻辑阵列维数 *输出：rank号 */int get_index(int row, int col, int sp) return (row+sp)%sp)*sp + (col+sp)%sp;/* *函数名：random_A_B *功能

17、：随机生成矩阵A和B */void random_A_B() int i,j; srand(unsigned int)time(NULL); /*设随机数种子*/*随机生成A,B,并初始化C*/ for(i=0; i<dg ; i+) for(j=0; j<dg ; j+) Aij = rand(); Bij = rand(); Cij = 0.0; /* 函数名：scatter_A_B * 功能：rank为0的处理器向其他处理器发送A、B矩阵的相关块 */void scatter_A_B() int i,j,k,l; int p_imin,p_imax,p_jmin,p_jmax

18、; for(k=0; k<p; k+) /*计算相应处理器所分得的矩阵块在总矩阵中的坐标范围*/ p_jmin = (k % sp ) * dl; p_jmax = (k % sp + 1) * dl-1; p_imin = (k - (k % sp)/sp * dl; p_imax = (k - (k % sp)/sp +1) *dl -1; l = 0; /*rank=0的处理器将A,B中的相应块拷至tmp_a,tmp_b，准备向其他处理器发送*/ for(i=p_imin; i<=p_imax; i+) for(j=p_jmin; j<=p_jmax; j+) tmp_

19、al = Aij; tmp_bl = Bij; l+; /*rank=0的处理器直接将自己对应的矩阵块从tmp_a,tmp_b拷至a,b*/ if(k=0) memcpy(a, tmp_a, dl2 * sizeof(float); memcpy(b, tmp_b, dl2 * sizeof(float); else /*rank=0的处理器向其他处理器发送tmp_a,tmp_b中相关的矩阵块*/ MPI_Send(tmp_a, dl2, MPI_FLOAT, k, 1, MPI_COMM_WORLD); MPI_Send(tmp_b, dl2, MPI_FLOAT, k, 2, MPI_CO

20、MM_WORLD); /* *函数名:init_alignment *功能:矩阵A和B初始对准 */void init_alignment() /*将A中坐标为(i,j)的分块A(i,j)向左循环移动i步*/ MPI_Sendrecv(a, dl2, MPI_FLOAT, get_index(my_row,my_col-my_row,sp), 1, tmp_a, dl2, MPI_FLOAT, get_index(my_row,my_col+my_row,sp), 1, MPI_COMM_WORLD, &status); memcpy(a, tmp_a, dl2 * sizeof(fl

21、oat) ); /*将B中坐标为(i,j)的分块B(i,j)向上循环移动j步*/ MPI_Sendrecv(b, dl2, MPI_FLOAT, get_index(my_row-my_col,my_col,sp), 1, tmp_b, dl2, MPI_FLOAT, get_index(my_row+my_col,my_col,sp), 1, MPI_COMM_WORLD, &status); memcpy(b, tmp_b, dl2 * sizeof(float) );/* *函数名：main_shift *功能：分块矩阵左移和上移，并计算分块c */void main_shift

22、() int i,j,k,l; for(l=0; l<sp; l+) /*矩阵块相乘，c+=a*b */ for(i=0; i<dl; i+) for(j=0; j<dl; j+) for(k=0; k<dl; k+) ci*dl+j += ai*dl+k*bk*dl+j; /* 将分块a左移1位 */ MPI_Send(a , dl2, MPI_FLOAT, get_index(my_row, my_col-1, sp), 1, MPI_COMM_WORLD); MPI_Recv(a , dl2, MPI_FLOAT, get_index(my_row, my_col

23、+1, sp), 1, MPI_COMM_WORLD, &status); /* 将分块b上移1位 */ MPI_Send(b , dl2, MPI_FLOAT, get_index(my_row-1, my_col, sp), 1, MPI_COMM_WORLD); MPI_Recv(b , dl2, MPI_FLOAT, get_index(my_row+1, my_col, sp), 1, MPI_COMM_WORLD, &status); /* *函数名：collect_c *功能：rank为0的处理器从其余处理器收集分块矩阵c */void collect_C() i

24、nt i,j,i2,j2,k; int p_imin,p_imax,p_jmin,p_jmax; /* 分块矩阵在总矩阵中顶点边界值 */ /* 将rank为0的处理器中分块矩阵c结果赋给总矩阵C对应位置 */ for (i=0;i<dl;i+) for(j=0;j<dl;j+) Cij=ci*dl+j; for (k=1;k<p;k+) /*将rank为0的处理器从其他处理器接收相应的分块c*/ MPI_Recv(c, dl2, MPI_FLOAT, k, 1, MPI_COMM_WORLD, &status); p_jmin = (k % sp ) *dl; p_

25、jmax = (k % sp + 1) *dl-1; p_imin = (k - (k % sp)/sp *dl; p_imax = (k - (k % sp)/sp +1) *dl -1; i2=0; /*将接收到的c拷至C中的相应位置,从而构造出C*/ for(i=p_imin; i<=p_imax; i+) j2=0; for(j=p_jmin; j<=p_jmax; j+) Cij=ci2*dl+j2; j2+; i2+; /*函数名：print *功能：打印矩阵 *输入：指向矩阵指针的指针，字符串 */void print(float *m,char *str) int

26、i,j; printf("%s",str); /*打印矩阵m*/ for(i=0;i<dg;i+) for(j=0;j<dg;j+) printf("%15.0f ",mij); printf("n"); printf("n");/* *函数名：main *功能：主过程，Cannon算法，矩阵相乘 *输入：argc为命令行参数个数，argv为每个命令行参数组成的字符串数组 */int main(int argc, char *argv) int i; MPI_Init(&argc, &a

27、rgv); /* 启动MPI计算 */ MPI_Comm_size(MPI_COMM_WORLD, &p); /* 确定处理器个数 */ MPI_Comm_rank(MPI_COMM_WORLD, &my_rank); /* 确定各自的处理器标识符 */ sp = sqrt(p); /* 确保处理器个数是完全平方数，否则打印错误信息，程序退出 */ if (sp*sp != p) if (my_rank = 0) printf("Number of processors is not a quadratic number!n"); MPI_Finalize(

28、); exit(1); if (argc != 2) if (my_rank = 0) printf("usage: mpirun -np ProcNum cannon MatrixDimensionn"); MPI_Finalize(); exit(1); dg = atoi(argv1); /* 总矩阵维数 */ dl = dg / sp; /* 计算分块矩阵维数 */ dl2 = dl * dl; /* 计算处理器在逻辑阵列中的坐标 */ my_col = my_rank % sp ; my_row = (my_rank-my_col) / sp ; /* 为a、b、

29、c分配空间 */ a = (float *)malloc( dl2 * sizeof(float) ); b = (float *)malloc( dl2 * sizeof(float) ); c = (float *)malloc( dl2 * sizeof(float) ); /* 初始化c */ for(i=0; i<dl2 ; i+) ci = 0.0; /* 为tmp_a、tmp_b分配空间 */ tmp_a = (float *)malloc( dl2 * sizeof(float) ); tmp_b = (float *)malloc( dl2 * sizeof(float

30、) ); if (my_rank = 0) /* rank为0的处理器为A、B、C分配空间 */ starttime = MPI_Wtime(); A = (float *)malloc( dg * sizeof(float*) ); B = (float *)malloc( dg * sizeof(float*) ); C = (float *)malloc( dg * sizeof(float*) ); for(i=0; i<dg; i+) Ai = (float *)malloc( dg * sizeof(float) ); Bi = (float *)malloc( dg * s

31、izeof(float) ); Ci = (float *)malloc( dg * sizeof(float) ); random_A_B(); /* rank为0的处理器随机化生成A、B矩阵 */ scatter_A_B(); /* rank为0的处理器向其他处理器发送A、B矩阵的相关块 */ else /* rank不为0的处理器接收来自rank为0的处理器的相应矩阵分块 */ MPI_Recv(a, dl2, MPI_FLOAT, 0 , 1, MPI_COMM_WORLD, &status); MPI_Recv(b, dl2, MPI_FLOAT, 0 , 2, MPI_CO

32、MM_WORLD, &status); init_alignment(); /* A、B矩阵的初始对准 */ main_shift(); /* 分块矩阵左移、上移, cannon算法的主过程 */ if(my_rank = 0) collect_C(); /* rank为0的处理器从其余处理器收集分块矩阵c */ print(A,"random matrix A : n"); /* 打印矩阵A */ print(B,"random matrix B : n"); /* 打印矩阵B */ print(C,"Matrix C = A * B

33、 : n"); /* 打印矩阵C */ else MPI_Send(c,dl2,MPI_FLOAT,0,1,MPI_COMM_WORLD); /* rank不为0的处理器向rank为0的处理器发送矩阵块c */ MPI_Barrier(MPI_COMM_WORLD); /* 同步所有处理器 */if(my_rank = 0) time1 = MPI_Wtime();printf("%s",time1-starttime); /* 运行时间*/ MPI_Finalize(); /* 结束MPI计算 */ return 0;四、结果分析和总结由第三部分数据结果可以看出，多个处理器并行算法比单个处理器串行程序在计算大量数据时要快得多。分析如下：设，相乘，结果矩阵：，其中，记一次乘法和加法为一个运算