MPI用于矩阵乘积示例

1.准备

使用MPI做并行计算时，根据程序的具体要求，可按任务进行分配或数据进行分配。根据矩阵乘积的特点，这里按数据进行分配，即每个计算机节点计算不同的数据，由于矩阵数据的特点，这里按行进行数据分块。因为本人使用的是C语言，数组在C语言的表示下行间数据地址是连续的(注：若是Fortran语言，则列是连续的)。

2.mpi程序的框架

mpi程序运行是靠输入dos命令执行的，因此，mpi程序一般都在main函数内，也即程序入口函数中。一般都会有MPI_Init、MPI_Comm_rank、MPI_Comm_size、MPI_Finalize这四个函数。用法如下：

记得要加上mpi.lib

 1 #include “mpi.h”
 2 
 3 int main(int argc,char *argv[])
 4 {
 5     int myid,numprocs;
 6     MPI_Init(&argc,&argv);//MPI Initialize
 7     MPI_Comm_rank(MPI_COMM_WORLD,&myid);//获得当前进程号
 8     MPI_Comm_size(MPI_COMM_WORLD,&numprocs);//获得进程个数
 9     //mpi计算过程
10     MPI_Finalize();//结束
11 }

3.矩阵乘法

矩阵乘法在于对矩阵进行分块，然后交由各进程执行，最后将计算结果传递给主进程。

假设是M*N，计算前，将矩阵N发送给所有从进程，然后将矩阵M分块，将M中数据按行分给各从进程，在从进程中计算M中部分行数据和N的乘积，最后将结果发送给主进程。这里为了方便，有多少进程，就将M分了多少块，除最后一块外的其他数据块大小都相等，最后一块是剩下的数据，大小大于等于其他数据块大小，因为矩阵行数不一定整除进程数。最后一块数据在主进程中计算，其他的在从进程中计算。

定义两个矩阵M和N，N所有进程都需要，M可以只在主进程中定义。其他的变量视主进程和从进程需要按要求定义在合适的位置。

代码如下，包括矩阵初始化，数据传递，矩阵乘积计算等。

  1 void matgen(float* a,int Width);
  2 //产生矩阵
  3 void matgen(float* a,int Width)
  4 {
  5     int i,j;
  6     for (i=0;i<Width;i++)
  7     {
  8         for (j=0;j<Width;j++)
  9         {
 10             //a[i*Width+j]=(float)rand()/RAND_MAX + (float)rand()/(RAND_MAX*RAND_MAX); //产生矩阵,矩阵中元素0~1
 11             a[i*Width+j]=1.00;
 12         }
 13     }
 14 }
 15 
 16 void main(int argc,char *argv[])
 17 {
 18     float *M,*N,*P,*buffer,*ans;
 19     int Width=1000;
 20     int myid,numprocs;
 21     MPI_Status status;
 22 
 23     MPI_Init(&argc,&argv);//MPI Initialize
 24     MPI_Comm_rank(MPI_COMM_WORLD,&myid);//获得当前进程号
 25     MPI_Comm_size(MPI_COMM_WORLD,&numprocs);//获得进程个数
 26 
 27     int line = Width/numprocs;//将数据分为(进程数)个块,主进程也要处理数据
 28     M = (float*)malloc(sizeof(float)*Width*Width);
 29     N = (float*)malloc(sizeof(float)*Width*Width);
 30     P = (float*)malloc(sizeof(float)*Width*Width);
 31         //缓存大小大于等于要处理的数据大小，大于时只需关注实际数据那部分
 32     buffer = (float*)malloc(sizeof(float)*Width*line);//数据分组大小
 33     ans = (float*)malloc(sizeof(float)*Width*line);//保存数据块结算的结果
 34 
 35     //主进程对矩阵赋初值，并将矩阵N广播到各进程,将矩阵M分组广播到各进程
 36     if (myid==0)
 37     {
 38         //矩阵赋初值
 39         matgen(M,Width);
 40         matgen(N,Width);
 41         //将矩阵N发送给其他从进程
 42         //MPI_Bcast(N,Width*Width,MPI_FLOAT,0,MPI_COMM_WORLD);
 43         for (int i=1;i<numprocs;i++)
 44         {
 45             MPI_Send(N,Width*Width,MPI_FLOAT,i,0,MPI_COMM_WORLD);
 46         }
 47         //依次将M的各行发送给各从进程
 48         for (int m=1;m<numprocs;m++)
 49         {
 50             /*
 51             //按行分组，得到buffer，发送给从进程
 52             for (int i=(m-1)*line;i<m*line;i++)
 53             {
 54                 for (int j=0;j<Width;j++)
 55                 {
 56                     buffer[(i-(m-1)*line)*Width+j]=M[i*Width+j];
 57                 }            
 58             }    
 59             //经数据buffer发送到进程m
 60             MPI_Send(buffer,line*Width,MPI_FLOAT,m,m,MPI_COMM_WORLD);
 61             */
 62             MPI_Send(M+(m-1)*line*Width,Width*line,MPI_FLOAT,m,1,MPI_COMM_WORLD);
 63         }   
 64         //接收从进程计算的结果
 65         for (int k=1;k<numprocs;k++)
 66         {
 67             MPI_Recv(ans,line*Width,MPI_FLOAT,k,3,MPI_COMM_WORLD,&status);
 68             //将结果传递给数组P
 69             for (int i=0;i<line;i++)
 70             {
 71                 for (int j=0;j<Width;j++)
 72                 {
 73                     P[((k-1)*line+i)*Width+j]=ans[i*Width+j];
 74                 }
 75             }
 76         }
 77         //计算M剩下的数据
 78         for (int i=(numprocs-1)*line;i<Width;i++)
 79         {
 80             for (int j=0;j<Width;j++)
 81             {
 82                 float temp=0.0;
 83                 for (int k=0;k<Width;k++)
 84                     temp += M[i*Width+k]*N[k*Width+j];
 85                 P[i*Width+j]=temp;
 86             }
 87         }
 88         //测试结果,计算一行的和
 89         float sum1=0;
 90         float sum2=0;
 91         for (int i=0;i<Width;i++)
 92         {
 93             sum1 += M[i];
 94             sum2 += P[600*Width+i];
 95         }
 96         printf("sum1=%f  sum2=%f\n",sum1,sum2);
 97         //统计时间
 98         double clockEnd = (double)clock();
 99         printf("myid:%d time:%.2fs\n",myid,(clockEnd-clockStart)/1000);//结果测试
100     }
101 
102     //其他进程接收数据，计算结果后，发送给主进程
103     else
104     {
105         //接收广播的数据(矩阵N)
106         //MPI_Bcast(N,Width*Width,MPI_FLOAT,0,MPI_COMM_WORLD);
107         MPI_Recv(N,Width*Width,MPI_FLOAT,0,0,MPI_COMM_WORLD,&status);
108         
109         MPI_Recv(buffer,Width*line,MPI_FLOAT,0,1,MPI_COMM_WORLD,&status);
110         //计算乘积结果，并将结果发送给主进程
111         for (int i=0;i<line;i++)
112         {
113             for (int j=0;j<Width;j++)
114             {
115                 float temp=0.0;
116                 for(int k=0;k<Width;k++)
117                     temp += buffer[i*Width+k]*N[k*Width+j];
118                 ans[i*Width+j]=temp;
119             }
120         }
121         //将计算结果传送给主进程
122         MPI_Send(ans,line*Width,MPI_FLOAT,0,3,MPI_COMM_WORLD);
123     }
124 
125     MPI_Finalize();//结束
126 
127     return 0;
128 }

可以对上述结果进行测试，测试时为了方便可以不用每次都在dos终端下输入命令，建立一个.bat的批处理命令即可，方法如下：

新建一个文本文件，在里面加上执行命令，诸如：mpiexec -n 5 XXX.exe，然后将文件重命名为xxx.bat，双击该bat文件即可执行mpi程序，但是会看不到输出结果，此方法在带有GUI的程序时很有用，运行其他程序或要对运行参数进行修改时，只需要编辑bat文件即可。相对于在命令提示符下可谓方便许多。

对于上面的程序其实还有优化的地方，仔细分析，可知，进程中需要的数据如下：

主进程：M、N、P(存储结果)、line

从进程：N、buffer(存储M分块的数据)、ans(存储计算结果)、line(计算M的行数)

buffer直接接受M分块的数据，主进程中矩阵P直接接受从进程ans传递的结果值。

如下：

1 for (int k=1;k<numprocs;k++)
2 {                    
      MPI_Recv(P+(k-1)*line*Width,line*Width,MPI_FLOAT,k,3,MPI_COMM_WORLD,&status);
3 }

4.题外

除了该方法外，还有一种就是计算都在从进程中计算，主进程只负责数据的管理。这样计算所花的时间要比上面方法的时间短，具体请自行验证。