【MPI学习5】MPI并行程序设计模式：组通信MPI程序设计
相关章节：第13章组通信MPI程序设计。
MPI组通信与点到点通信的一个重要区别就是：组通信需要特定组内所有成员参与，而点对点通信只涉及到发送方和接收方。
由于需要组内所有成员参与，因此也是一种比较复杂的通信方式。程序员在设计组通信语句的时候，需要同时考虑两点：
a. 程序运行起来之后，当前正在运行的进程的行为方式
b. 将组通信作为一个整体，考虑所有进程的行为方式
（1）概述
组通信从功能上实现了三个方面：
a. 通信：完成组内数据传输（广播、收集、散发、组收集、全互换各种数据交换传输方式）
b. 同步：能够让组内所有进程在特定的位置上执行进度上取得一致，组通信虽然可以有同步功能，但是并不意味这同步一定发生（MPI_Barrier同步函数）
c. 计算：通过给定数据的OP归约操作完成（分为MPI预定义OP归约操作、用户自定义OP归约操作）
下面分别阐述上述三个方面的相关代码示例，采用的方法是把书上的代码段扩充成完整的可执行的程序来体会。
（2）各部分代码
下面阐述几个代码示例，把通信、同步、计算的例子都融入到具体的代码中。
示例1（Bcast广播通信）
从root进程读入一个数，并Bcast广播给组内所有进程并打印到终端上，代码如下：
#include "mpi.h"
#include <stdio.h>
#include <stdlib.h>
int main(int argc, char *argv[])
{
int rank, value=0;
MPI_Init(&argc, &argv);
MPI_Comm_rank(MPI_COMM_WORLD, &rank);
while(value>=0)
{
if (0==rank) {
scanf("%d", &value);
}
MPI_Bcast(&value, 1, MPI_INT, 0, MPI_COMM_WORLD);
printf("Process %d got %d\n", rank, value);
}
MPI_Finalize();
return 0;
}
代码执行结果如下：
上述代码中涉及到了Bcast函数的基本用法，在每个进程中都要出现，然后通过函数参数中的root进程号，进而区分哪个进程是发送方。
示例2（Alltoall全互换通信）
全互换。全互换是形式上最复杂的一种组通信方式，内容上是各种通信形式的全集。
全互换通信模式中，在组内所有通信涉及到的进程中：
a. 每个进程发送给其他进程的数据是不同的
b. 每个进程从其他进程获得数据也是不同的
这样描述还是太抽象，通过具体的图来解释一下：
上图描述了全互换发生前和发生后，数据变化的情况：
纵轴代表不同进程编号；左图的横轴代表发送缓冲区的段位大小，右图的横轴代表接受缓冲区的大小。
a. 第i个进程的发送缓冲区中的第j段数据（即左图Aij）代表的意义是：第i个进程发送给第j个进程的数据
b. 第j个进程的接受缓冲区中的第i段数据（即右图的Aji）代表的意义是：第j个进程从第i个进程中接收到的数据
结合a、b两点，我们可以知道全互换这种方式类似：将“进程-发送缓冲区”矩阵，转置成“进程-接受缓冲区”矩阵
对于同一对（i,j），左图的Aij的大小要求一定与右图的Aji的大小一样。
比如第0个进程发送给第2个进程的数据（即A02)，必须与第2个进程给第0个进程留出来的接受缓冲区大小一致（即A20），满足这样的条件，就实现了全互换。
每个进程在全通信中既向所有进程发送数据，又同时从所有进程接受数据，但是需要注意上述的发送、接受缓冲区大小匹配。
看如下的代码：
1 #include "mpi.h"
2 #include <stdlib.h>
3 #include <stdio.h>
4 #include <string.h>
5 #include <errno.h>
6
7 int main(int argc, char *argv[])
8 {
9 int rank, size;
10 int chunk = 2; // 发送到一个进程的数据块大小
11 int i,j;
12 int *sbuf, *rbuf;
13
14 MPI_Init(&argc, &argv);
15 MPI_Comm_rank(MPI_COMM_WORLD, &rank);
16 MPI_Comm_size(MPI_COMM_WORLD, &size);
17 // 申请发送缓冲区 和 接收缓冲区
18 sbuf = (int*)malloc(size*chunk*sizeof(int)); // 申请发送缓冲区
19 if (!sbuf) {
20 perror("can't allocate send buffer");
21 MPI_Abort(MPI_COMM_WORLD, EXIT_FAILURE);
22 }
23 rbuf = (int*)malloc(size*chunk*sizeof(int)); // 申请接收缓冲区
24 if (!rbuf) {
25 perror("can't allocate recv buffer");
26 free(sbuf);
27 MPI_Abort(MPI_COMM_WORLD, EXIT_FAILURE);
28 }
29 // 设置每个进程的发送缓冲区 和 接收缓冲区的数据
30 for(i=0; i<size; i++) // i代表进程编号
31 {
32 for(j=0; j<chunk; j++) // j代表myid发送给进程i的第j个数据
33 {
34 sbuf[i*chunk+j] = rank+i*chunk+j; // 设置发送缓冲区的数据
35 printf("myid = %d, send to id = %d, data[%d] = %d\n",rank ,i ,j, sbuf[i*chunk+j]);
36 rbuf[i*chunk+j] = 0; // 接收缓冲区清零
37 }
38 }
39 // 执行all to all调用
40 MPI_Alltoall(sbuf, chunk, MPI_INT, rbuf, chunk, MPI_INT, MPI_COMM_WORLD);
41 // 打印接收缓冲区从其他进程接收的数据
42 for(i=0; i<size; i++)
43 {
44 for(j=0; j<chunk; j++)
45 {
46 printf("myid = %d, recv from id = %d, data[%d] = %d\n",rank ,i ,j, rbuf[i*chunk+j]);
47 }
48 }
49 free(rbuf);
50 free(sbuf);
51 MPI_Finalize();
52 return 0;
53 }
程序执行结果如下：
一共有2个进程，每个进程共有4个数据；可以看到执行全互换后的效果。
示例3
利用近似积分公式求圆周率π的大小。（Bcast广播、Barrier同步、REDUCE归约）
圆周率的大小可以用一个积分表达式来求解，求这个积分的过程可以用MPI组通信技术实现
代码如下：
1 #include "mpi.h"
2 #include <stdio.h>
3 #include <stdlib.h>
4
5 double f(double);
6
7 int main(int argc, char *argv[])
8 {
9 int rank, size;
10 int n; // 细分出来的小矩阵个数
11 double pi; // 存放pi的计算值
12
13 MPI_Init(&argc, &argv);
14 MPI_Comm_size(MPI_COMM_WORLD, &size);
15 MPI_Comm_rank(MPI_COMM_WORLD, &rank);
16 if (0==rank) {
17 printf("number of small rectangle:\n");
18 scanf("%d",&n);
19 }
20 MPI_Bcast(&n, 1, MPI_INT, 0, MPI_COMM_WORLD); // 将n值广播到各个进程中
21 int i;
22 double x, h, sum=0.0;
23 for (i=rank+1; i<=n; i+=size)
24 {
25 x = (i-0.5)/(double)n;
26 sum += f(x);
27 }
28 sum = sum/n;
29 MPI_Reduce(&sum, &pi, 1, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD);
30 if (0==rank) {
31 printf("Approximately of pi is : %.16f\n",pi);
32 fflush(stdout);
33 }
34 MPI_Finalize();
35 }
36
37 double f(double x) { return 4.0/(1.0+x*x); }
上述代码是按照如下思路实现的：
比如总共要把积分区间按X轴划分成1000份，即计算任务是求1000个小矩形的面积，
而MPI计算进程共有5个；则为了负载均衡，考虑让每个计算节点平均搞定200个小矩形的面积。
代码技巧：
这里有个需要处理的地方是：如何告诉每个计算节点（进程），需要计算哪些小矩形的面积呢？
一种直观的思路是，让第0个进程处理0~199小矩形，让第1个进程处理200~399小矩形...，以此类推。
上述的思路没有错误，但是对边界条件的处理可能麻烦一些：因为用这种连续区间的任务分配方式，就需要知道提前知道头和尾各在哪里，需要额外处理。
书上给的代码技巧是通过间隔划分的方式，比如有5个计算进程，第0个进程计算0，5，10，...这些矩形面积，第1个进程计算1，6，11,...这些矩形面积。
这样一来，任务分配的代码就简单了，一个循环直接搞定，而且不同进程不用区别对待。
上述代码执行结果如下：
小矩形个数越多，计算的值约逼近真实值。
示例4
有一个大数据集，每个计算进程中各有数据集中的一部分，现在想让每个进程都有全部的数据集（Bcast广播、Barrier同步）
代码如下：
1 #include "mpi.h"
2 #include <stdlib.h>
3 #include <stdio.h>
4
5 int main(int argc, char *argv[])
6 {
7 int rank,size,i;
8 int *table;
9 int errors = 0;
10 MPI_Aint address;
11 MPI_Datatype type, newtype;
12 int lens;
13
14 MPI_Init(&argc, &argv);
15 MPI_Comm_rank(MPI_COMM_WORLD, &rank);
16 MPI_Comm_size(MPI_COMM_WORLD, &size);
17 table = (int*)calloc(size, sizeof(int));
18 table[rank] = rank + 1; // 准备要广播的数据
19 MPI_Barrier(MPI_COMM_WORLD);
20 for(i=0; i<size; i++) MPI_Bcast(&table[i], 1, MPI_INT, i, MPI_COMM_WORLD); // 每个进程只当一次root
21 for(i=0; i<size; i++)
22 {
23 if (table[i]!=i+1) {
24 errors++;
25 }
26 }
27 MPI_Barrier(MPI_COMM_WORLD);
28 if (0==rank) {
29 for(i=0;i<size;i++) printf("table[%d]:%d\n",i,table[i]);
30 }
31 MPI_Finalize();
32 }
代码执行结果如下：
代码分析：
a. line19执行了第一次Barrier同步，目的是让各个进程都准备好各自拥有数据集的一部分（这里同步的目的是，让各个进程一方面table分配好空间，并且准备好数据）
b. line27执行了第二次Barrier同步，目的是让Bcast都完成，然后在root进程中查看table数据是否都发送过去了
其中line20的代码技巧也值得学习，每个进程只当一次root，而其余的时候全都是当配角，通信中的接收一方
示例5
组通信中的死锁。组通信的相同函数有可能是同步的，也有可能是异步的，这个受软硬件平台影响；为了构建可移植性好的代码，在设计组通信程序的时候，无论是异步还是同步的，都应该保证程序不出现死锁。下面用几个例子解释死锁出现的情况。
代码1. 同一个通信域中的死锁示例
1 #include "mpi.h"
2 #include <stdio.h>
3 #include <stdlib.h>
4 #include <unistd.h>
5
6 #define LEN 10
7
8 int main(int argc, char *argv[])
9 {
10 int rank,size;
11
12 MPI_Init(&argc, &argv);
13 MPI_Comm_rank(MPI_COMM_WORLD, &rank);
14 MPI_Comm_size(MPI_COMM_WORLD, &size);
15
16 if(2!=size) MPI_Abort(MPI_COMM_WORLD, EXIT_FAILURE);
17 int* buf1 = (int*)malloc(sizeof(int)*LEN);
18 int* buf2 = (int*)malloc(sizeof(int)*LEN);
19 buf1[0] = 1;
20 buf2[0] = 1;
21 if (0==rank) {
22 MPI_Bcast(buf1,LEN,MPI_INT,0,MPI_COMM_WORLD);
23 MPI_Bcast(buf2,LEN,MPI_INT,1,MPI_COMM_WORLD);
24 printf("0 done\n");
25 }
26 if (1==rank) {
27 MPI_Bcast(buf2,LEN,MPI_INT,1,MPI_COMM_WORLD);
28 MPI_Bcast(buf1,LEN,MPI_INT,0,MPI_COMM_WORLD);
29 printf("1 done\n");
30 }
31 free(buf1);
32 free(buf2);
33 MPI_Finalize();
34 }
代码执行结果如下：
可以看到并没有出现deadlock。原因是由于LEN定义的发送和接收数据长度太小，MPI给Bcast执行方式选择的是缓存非同步方式，即Bcast不必等到所有需要参加组通信的进程都完成了再返回。
如果我们把LEN定义的长度改为100000，则执行结果如下：
执行结果表明deadlock出现了。原因是由于LEN定义的发送和接受的数据长度较大，MPI给Bcast选择的执行方式是同步方式，即Bcast必须等到所有需要参加广播通信的进程都完成了，才能够正确返回。
书上并没有提Bcast是异步还是同步的事情，我在stackoverflow上提问才知道了上面的道理（http://stackoverflow.com/questions/35628524/why-this-mpi-bcast-related-code-not-deadlock）。感谢stackoverflow上的认真回答。
代码2. 不同通信域中的死锁示例
这部分内容与后面第15章的进程组和通信域知识相关，需要调换顺序提前补充。
总共有3个进程，三个进程在MPI_COMM_WORLD中的rank分别是{0，1，2}
现在构造三个新的通信域comm0 comm1 comm2
comm0中包含MPI_COMM_WORLD中rank为0，1的进程
comm1中包含MPI_COMM_WORLD中rank为1，2的进程
comm2中包含MPI_COMM_WORLD中rank为2，0的进程
需要注意是，虽然进程客观上只有3个，但是相同的进程在不同的通信域中的rank号是不同的。
下面的代码构造了一个循环依赖的deadlock
1 #include "mpi.h"
2 #include <stdio.h>
3 #include <stdlib.h>
4
5 #define LEN 100000
6
7 int main(int argc, char *argv[])
8 {
9 MPI_Init(&argc, &argv);
10 int rank, size;
11 MPI_Comm_size(MPI_COMM_WORLD, &size);
12 if (3!=size) {
13 printf("wrong proc num\n");
14 MPI_Abort(MPI_COMM_WORLD, EXIT_FAILURE);
15 }
16 const int ranks[3] = {0,1,2};
17 MPI_Group world_group;
18 MPI_Comm_group(MPI_COMM_WORLD, &world_group);
19 // 构造三个新的进程组对象
20 MPI_Group group0, group1, group2;
21 const int ranks0[2] = {0,1};
22 const int ranks1[2] = {1,2};
23 const int ranks2[2] = {2,0};
24 MPI_Group_incl(world_group, 2, ranks0, &group0);
25 MPI_Group_incl(world_group, 2, ranks1, &group1);
26 MPI_Group_incl(world_group, 2, ranks2, &group2);
27
28 // 利用三个进程对象构造通信域
29 MPI_Comm comm0, comm1, comm2;
30 MPI_Comm_create(MPI_COMM_WORLD, group0, &comm0);
31 MPI_Comm_create(MPI_COMM_WORLD, group1, &comm1);
32 MPI_Comm_create(MPI_COMM_WORLD, group2, &comm2);
33 MPI_Comm_rank(MPI_COMM_WORLD, &rank);
34 int *buf1 = (int*)malloc(sizeof(int)*LEN);
35 int *buf2 = (int*)malloc(sizeof(int)*LEN);
36 buf1[0] = 1;
37 buf2[0] = 1;
38 int sub_rank;
39 if (0==rank) {
40 MPI_Group_translate_ranks(world_group,1,&ranks[0],group0,&sub_rank);
41 MPI_Bcast(buf1, LEN, MPI_INT, sub_rank, comm0);
42 MPI_Group_translate_ranks(world_group,1,&ranks[2],group2,&sub_rank);
43 MPI_Bcast(buf2, LEN, MPI_INT, sub_rank, comm2);
44 }
45 if (1==rank) {
46 MPI_Group_translate_ranks(world_group,1,&ranks[1],group1,&sub_rank);
47 MPI_Bcast(buf1, LEN, MPI_INT, sub_rank, comm1);
48 MPI_Group_translate_ranks(world_group,1,&ranks[0],group0,&sub_rank);
49 MPI_Bcast(buf2, LEN, MPI_INT, sub_rank, comm0);
50 }
51 if (2==rank) {
52 MPI_Group_translate_ranks(world_group,1,&ranks[2],group2,&sub_rank);
53 MPI_Bcast(buf1, LEN, MPI_INT, sub_rank, comm2);
54 MPI_Group_translate_ranks(world_group,1,&ranks[1],group1,&sub_rank);
55 MPI_Bcast(buf2, LEN, MPI_INT, sub_rank, comm1);
56 }
57 MPI_Finalize();
58 }
代码执行结果是deadlock，这个deadlock本身不难理解，就是comm0依赖comm1，comm1依赖comm2，comm2依赖comm0，循环依赖造成了死锁，互相都等着。
书上只给出来了关键代码的原型，并没有给出来通信域是怎么构造的，Bcast中的sub_rank是怎么得来的，需要补充一些进程组和通信域的知识。
进程组和通信域使得MPI程序与传统的串行程序的设计思路上有很大不同，初学甚至有些别扭，我也没有全部领会。
额外参考了https://www.rc.usf.edu/tutorials/classes/tutorial/mpi/chapter9.html
只能先把与这个例子相关的知识记录下来：
a. 通信域（Communicator）是大的概念，进程组是通信域的一个属性
b. 一个通信域对应一个进程组
c. 一个进程可以在多个进程组中，因此也可以在多个通信域中
d. 通过进程组来构造通信域是通信域的一种生成方式
e. 相同的进程在不同的通信域中的rank是不同的，可以通过函数来查阅这种对应关系
f. 比如，调用MPI_Comm_create函数，利用进程组{0，1}构造一个通信域comm0；那么进程2也会执行这条语句，则进程2中对应的comm0就是MPI_COMM_NULL
关于上面提到的f，我用http://mpitutorial.com/tutorials/introduction-to-groups-and-communicators/里面的一个代码示例进行了验证，加深理解：
1 #include "mpi.h"
2 #include <stdio.h>
3 #include <stdlib.h>
4
5 int main(int argc, char *argv[])
6 {
7 MPI_Init(&argc, &argv);
8 int world_rank, world_size;
9 MPI_Comm_rank(MPI_COMM_WORLD, &world_rank);
10 MPI_Comm_rank(MPI_COMM_WORLD, &world_size);
11
12 MPI_Group world_group;
13 MPI_Comm_group(MPI_COMM_WORLD, &world_group);
14
15 int n = 7;
16 const int ranks[7] = {1,2,3,5,7,11,13};
17
18 MPI_Group prime_group;
19 // 构造新的进程组 这个进程组是客观的 不管当前进程是哪个
20 MPI_Group_incl(world_group, 7, ranks, &prime_group);
21
22 MPI_Comm prime_comm;
23 // 如果当前进程不在prime_group中 则prime_com就是MPI_COMM_NULL
24 MPI_Comm_create_group(MPI_COMM_WORLD, prime_group, 0, &prime_comm);
25
26 int prime_rank = -1, prime_size = -1;
27 if (MPI_COMM_NULL!=prime_comm) {
28 MPI_Comm_rank(prime_comm, &prime_rank);
29 MPI_Comm_size(prime_comm, &prime_size);
30 }
31
32 printf("WORLD RANK/SIZE: %d/%d \t PRIME RANK/SIZE: %d/%d\n",world_rank, world_size, prime_rank, prime_size);
33 if(MPI_COMM_NULL!=prime_comm) MPI_Barrier(prime_comm);
34 MPI_Barrier(MPI_COMM_WORLD);
35 if(MPI_GROUP_NULL!=world_group) MPI_Group_free(&world_group);
36 if(MPI_GROUP_NULL!=prime_group) MPI_Group_free(&prime_group);
37 if(MPI_COMM_NULL!=prime_comm) MPI_Comm_free(&prime_comm);
38 MPI_Finalize();
39 }
代码执行结果如下：
从原来进程中扣出来几个进程构造新进程；对于没被扣出来的那些进程，执行到MPI_Comm_create这里的时候，prime_comm就是MPI_COMM_NULL。
示例6
归约操作。这部分代码都是书上的内容，认真看书即可。
代码1 MINLOC和MAXLOC
有时候不仅需要知道一堆数据中的最小（最大）值，而且需要知道最小最大值对应的编号（或进程号）是多少。
通过MPI提供的MPI_REDUCE归约函数，以及预定义的归约操作MPI_MINLOC活MPI_MAXLOC即可方便实现。
1 #include "mpi.h"
2 #include <stdio.h>
3 #include <stdlib.h>
4 #include <time.h>
5
6 int main(int argc, char *argv[])
7 {
8 MPI_Init(&argc, &argv);
9 int i, myrank, root;
10 root = 0;
11 double ain[3]; // 各个进程的输入数据
12 double aout[3]; // 只有root进程有用 用于存放规约后的数字的具体值
13 int ind[3]; // 只有root进程有用 用于存放规约后的数字的编号
14 struct{
15 double val;
16 int rank;
17 } in[3], out[3]; // 规约操作的发送缓冲区和接收缓冲区
18 MPI_Comm_rank(MPI_COMM_WORLD, &myrank);
19 srand((unsigned int)(myrank+time(NULL)));
20 for(i=0;i<3;i++) ain[i] = (1.0*rand())/RAND_MAX;
21 for(i=0; i<3; i++){
22 in[i].val = ain[i];
23 in[i].rank = myrank;
24 }
25 printf("myrank:%d %.3f,%.3f,%.3f\n",myrank, ain[0],ain[1],ain[2]);
26 MPI_Reduce(in, out, 3, MPI_DOUBLE_INT, MPI_MAXLOC, root, MPI_COMM_WORLD);
27 if (myrank==root) {
28 printf("reduce result:\n");
29 for(i=0; i<3; i++){
30 aout[i] = out[i].val;
31 ind[i] = out[i].rank;
32 printf("%.3f of rank %d\n", aout[i], ind[i]);
33 }
34 }
35 MPI_Finalize();
36 }
代码执行结果如下：
可以看到这种归约操作：
a. 要求每个进程中的输入数据的大小是一样的（在这个例子中都是长度为3的double数组），
b. 归约的对象是不同进程的数组中相同位置的数据
另外，也允许用户自定义归约操作，如下面代码定义了复数相乘的归约操作（书上只给了代码框架，具体细节中的问题参考http://stackoverflow.com/questions/9285442/mpi-get-processor-with-minimum-value）。
1 #include "mpi.h"
2 #include <stdio.h>
3 #include <stdlib.h>
4 #include <time.h>
5
6 typedef struct{
7 double real;
8 double imag;
9 }Complex;
10
11 void multiplyOP(Complex *, Complex *, int *, MPI_Datatype *);
12
13 int main(int argc, char *argv[])
14 {
15 MPI_Init(&argc, &argv);
16 int root = 0;
17 int rank;
18 MPI_Comm_rank(MPI_COMM_WORLD, &rank);
19 Complex input[2], output[2];
20 if (0==rank) {
21 input[0].real = 1; input[0].imag = 1;
22 input[1].real = 1; input[1].imag = 2;
23 }
24 if (1==rank) {
25 input[0].real = 1; input[0].imag = -1;
26 input[1].real = 1; input[1].imag = 2;
27 }
28 // 在MPI中注册和构造复数类型
29 MPI_Datatype ctype;
30 MPI_Type_contiguous(2, MPI_DOUBLE, &ctype); // 在MPI中构造复数结构体 连续的两个DOUBLE类型
31 MPI_Type_commit(&ctype); // 在MPI中注册刚构造的复数结构体
32 // 在MPI中构造自定义归约操作
33 MPI_Op myOp;
34 MPI_Op_create((MPI_User_function*)multiplyOP, 1, &myOp); // 生成用户定义的乘积操作
35 MPI_Reduce(input, output, 2, ctype, myOp, root, MPI_COMM_WORLD); // 执行复数乘积操作
36 // 在root中打印结果
37 if (root==rank) {
38 printf("reduce result of 0 is : %1.0f+(%1.0f)i\n",output[0].real, output[0].imag);
39 printf("reduce result of 1 is : %1.0f+(%1.0f)i\n",output[1].real, output[1].imag);
40 }
41 MPI_Finalize();
42 }
43
44 // 计算复数的乘积
45 void multiplyOP(Complex *in, Complex *inout, int *len, MPI_Datatype *datatype)
46 {
47 int i;
48 Complex c;
49 for(i=0; i<*len; i++){
50 c.real = inout->real*in->real - inout->imag*in->imag;
51 c.imag = inout->real*in->imag + inout->imag*in->real;
52 *inout = c; // 把计算结果存回到inout的位置
53 in++;
54 inout++;
55 }
56 }
代码的执行结果如下：
