使用MPI进行分布式内存编程（2）
MPI的英文全称为message passing interface,中文名为消息传递接口，他不是一种新的语言，而是一个可以被C，C++，Fortran程序调用的库。
预备知识
1.编译与执行
使用类似此形式进行编译 mpicc -g -Wall -o mpi_hello mpi_hello.c 进行编译，mpicc为C语言的包装脚本(wrapper script)而非编译器(compilier)。
执行的话，可以使用 mpiexec -n <number of processers> ./mpi_hello 来执行，在一些超算平台，如天河二号上，也可以使用yhrun等命令运行
2.MPI Init and MPI Finalize
MPI_Init是告知MPI系统进行所有必要的初始化设置，如为消息缓冲区分配消息，为进程分配进程号等。
MPI_Finallize告知MPI系统已使用完毕，将MPI占用的资源释放。
一般的MPI程序框架如下：
. . . #include <mpi.h> . . . int main(int argc, char argv[]) f . . . /* No MPI calls before this */ MPI Init(&argc, &argv); . . . MPI Finalize(); /* No MPI calls after this */ . . . return 0;
｝
3.通信子，MPI_Comm_size和MPI_Comm_rank
通信子:一组可以相互发送消息的进程集合。
MPI_Comm_size：可以得到通信子的进程数
MPI_Comm_rank：可以得到通信子的进程号
4.MPI_Send和MPI_Recv
MPI_Send语法结构如下：
int MPI_Send( void msg_buf_p /* in */, int msg_size /* in */, MPI Datatype_msg_type /* in */, int dest /* in */, int tag /* in */, MPI_Comm_communicator /* in */);
msg_buf_p ---->指向包含消息内容的内存块的指针
msg_size ---->发生的数据量
Datatype_msg_type --->由于C语言之中的类型，如int，char等不能作为参数传递给函数。所以要使用这个参数作为媒介。
dest--->指定了需要接收信息的进程号
tag --->一个非负int类型数，用于区分看上去完全一样的信息。
MPI_Comm_communicator ---->通信子。
MPI_Recv语法结构如下：
int MPI_Recv( void msg_buf_p /* out */, int buf_size /* in */, MPI_Datatype buf_type /* in */, int source /* in */, int tag /* in */, MPI_Comm communicator /* in */, MPI_Status status_p /* out */);
msg_buf_p--->指向内存块
buf_size --->指定了内存块之中要储存对象的数量
buf_type --->说明对象的类型
source ---> 指定接受的信息应该由哪个进程发送而来
tag --->一个非负int类型数，用于区分看上去完全一样的信息,需要与发送消息的参数tag相匹配
communicator --->通信子，需要与发送消息的参数communicator 相匹配
status_p ----> 大部分情况下不进行调用。
实例见https://www.cnblogs.com/caishunzhe/p/12935439.html
为了编写能够使用scanf的MPI程序，我们根据进程号来选取分钟。0号进程复制读取数据，并将数据发送给其他进程。程序如下
pro3.5.c
#include<stdio.h> #include<string.h> #include<mpi.h> double f(double x) { return x*x+x*x*x+5; } double Trap(double left_endpt,double right_endpt,int trap_count,double base_len) { double estimate,x; int i; estimate=(f(left_endpt)+f(right_endpt))/2.0; //梯形面积 for(i =1;i<=trap_count-1;++i) { x=left_endpt+i*base_len; estimate+=f(x); } estimate=estimate*base_len; return estimate; } void Get_input( int my_rank, int comm_sz, double* a_p, double* b_p, int* n_p ) { int dest; if(my_rank==0) { printf("Enter a,b,and n\n"); scanf("%lf %lf %d",a_p,b_p,n_p); for(dest=1;dest<comm_sz;++dest) { MPI_Send(a_p,1,MPI_DOUBLE,dest,0,MPI_COMM_WORLD); MPI_Send(b_p,1,MPI_DOUBLE,dest,0,MPI_COMM_WORLD); MPI_Send(n_p,1,MPI_INT,dest,0,MPI_COMM_WORLD); } } else //rank!=0; { MPI_Recv(a_p,1,MPI_DOUBLE,0,0,MPI_COMM_WORLD,MPI_STATUS_IGNORE); MPI_Recv(b_p,1,MPI_DOUBLE,0,0,MPI_COMM_WORLD,MPI_STATUS_IGNORE); MPI_Recv(n_p,1,MPI_INT,0,0,MPI_COMM_WORLD,MPI_STATUS_IGNORE); } } int main() { /* int my_rank,comm_sz,n=1024,local_n; double a=0.0,b=3.0,h,local_a,local_b; */ int my_rank,comm_sz,n,local_n; double a,b,h,local_a,local_b; double local_int,total_int; int source; MPI_Init(NULL,NULL); //初始化 MPI_Comm_size(MPI_COMM_WORLD,&comm_sz); //返回进程数 MPI_Comm_rank(MPI_COMM_WORLD,&my_rank); //返回进程号 //新加入的函数，处理输入问题 Get_input(my_rank,comm_sz,&a,&b,&n); h=(b-a)/n; local_n=n/comm_sz; local_a=a+my_rank*local_n*h; local_b=local_a+local_n*h; local_int=Trap(local_a,local_b,local_n,h); if(my_rank!=0) { MPI_Send(&local_int,1,MPI_DOUBLE,0,0,MPI_COMM_WORLD); } else { total_int=local_int; for(source=1;source<comm_sz;source++) { MPI_Recv(&local_int,1,MPI_DOUBLE,source,0,MPI_COMM_WORLD,MPI_STATUS_IGNORE); //接受其他节点信息 total_int+=local_int; } } if(my_rank==0) { printf("With n=%d trapezoids,our estimated\n",n); printf("of the intergral from %f to %f=%.15e\n",a,b,total_int); } MPI_Finalize(); return 0; }
从0积分到50，划分的小间隔为100，1000,10000，理论上原函数为F(x)=1/3*x^3+1/4*x^4+5x ,积分值为F(50)-F(0)=1604416.6667
天河超算平台结果
可以看见，随着n的增大，划分的越来越细，结果越接近真实值。
实例上面是比较基本的实现，仔细分析之后发现其通信开销太大，0号进程需要负责所有进程的累加，如果又8个进程，那么就得加7次，随着进程的增长耗时线性增长
我们可以对其进行相应的改进，使用树型结构。
如下图的两种方法均可以，随进程数增长，时间为log（N）
我们选取第一种策略
#include<stdio.h> #include<string.h> #include<mpi.h> /*只考虑最简单的2的幂次方情况*/ double f(double x) { return x*x+x*x*x+5; } double Trap(double left_endpt,double right_endpt,int trap_count,double base_len) { double estimate,x; int i; estimate=(f(left_endpt)+f(right_endpt))/2.0; //梯形面积 for(i =1;i<=trap_count-1;++i) { x=left_endpt+i*base_len; estimate+=f(x); } estimate=estimate*base_len; return estimate; } void Get_input( int my_rank, int comm_sz, double* a_p, double* b_p, int* n_p ) { int dest; if(my_rank==0) { printf("Enter a,b,and n\n"); scanf("%lf %lf %d",a_p,b_p,n_p); for(dest=1;dest<comm_sz;++dest) { MPI_Send(a_p,1,MPI_DOUBLE,dest,0,MPI_COMM_WORLD); MPI_Send(b_p,1,MPI_DOUBLE,dest,0,MPI_COMM_WORLD); MPI_Send(n_p,1,MPI_INT,dest,0,MPI_COMM_WORLD); } } else //rank!=0; { MPI_Recv(a_p,1,MPI_DOUBLE,0,0,MPI_COMM_WORLD,MPI_STATUS_IGNORE); MPI_Recv(b_p,1,MPI_DOUBLE,0,0,MPI_COMM_WORLD,MPI_STATUS_IGNORE); MPI_Recv(n_p,1,MPI_INT,0,0,MPI_COMM_WORLD,MPI_STATUS_IGNORE); } } int main() { /* int my_rank,comm_sz,n=1024,local_n; double a=0.0,b=3.0,h,local_a,local_b; */ //comm_sz-->进程数 my_rank--->进程号 int my_rank,comm_sz,n,local_n; double a,b,h,local_a,local_b; double local_int,total_int; int source; int t; MPI_Init(NULL,NULL); //初始化 MPI_Comm_size(MPI_COMM_WORLD,&comm_sz); //返回进程数 MPI_Comm_rank(MPI_COMM_WORLD,&my_rank); //返回进程号 //新加入的函数，处理输入问题 Get_input(my_rank,comm_sz,&a,&b,&n); h=(b-a)/n; local_n=n/comm_sz; local_a=a+my_rank*local_n*h; local_b=local_a+local_n*h; local_int=Trap(local_a,local_b,local_n,h); /* if(my_rank!=0) { MPI_Send(&local_int,1,MPI_DOUBLE,0,0,MPI_COMM_WORLD); } else //my_rank=0 { total_int=local_int; for(source=1;source<comm_sz;source++) { MPI_Recv(&local_int,1,MPI_DOUBLE,source,0,MPI_COMM_WORLD,MPI_STATUS_IGNORE); //接受其他节点信息 total_int+=local_int; } } */ /* if(my_rank%2!=0)//奇数向小一位的偶数发送,然后退休 { MPI_Send(&local_int,1,MPI_DOUBLE,my_rank-1,0,MPI_COMM_WORLD); } else //偶数先接受奇数还可能有进一步的行动 { //接受奇数节点信息 total_int=local_int; MPI_Recv(&local_int,1,MPI_DOUBLE,my_rank+1,0,MPI_COMM_WORLD,MPI_STATUS_IGNORE); //接受其他节点信息 total_int+=local_int; local_int=total_int; for(t=2;t<comm_sz;t<<1) { for(source=t;source<comm_sz;source+=2*t) { if(my_rank==source) MPI_Send(&local_int,1,MPI_DOUBLE,my_rank-t,0,MPI_COMM_WORLD); } for(source=0;source<comm_sz;source+=2*t) { if(my_rank==source) { total_int=local_int; MPI_Recv(&local_int,1,MPI_DOUBLE,my_rank+t,0,MPI_COMM_WORLD,MPI_STATUS_IGNORE); //接受其他节点信息 total_int+=local_int; local_int=total_int; } } } } */ for(t=1;t<comm_sz;t*=2) { for(source=t;source<comm_sz;source+=2*t) { if(my_rank==source) { MPI_Send(&local_int,1,MPI_DOUBLE,my_rank-t,0,MPI_COMM_WORLD); break; } } for(source=0;source<comm_sz;source+=2*t) { if(my_rank==source) { total_int=local_int; MPI_Recv(&local_int,1,MPI_DOUBLE,my_rank+t,0,MPI_COMM_WORLD,MPI_STATUS_IGNORE); //接受其他节点信息 total_int+=local_int; local_int=total_int; break; } } } if(my_rank==0) { printf("With n=%d trapezoids,our estimated\n",n); printf("of the intergral from %f to %f=%.15e\n",a,b,total_int); //printf("comm_sz=%d",comm_sz); } MPI_Finalize(); return 0; }
天河超算平台结果
