MPI教程2

§点对点通信

MPI最基本的通信模式是在一对进程之间进行的消息收发操作：一个进程发送消息，另一个进程接收消息，这种通信方式称为点对点通信(point to point communications)。

MPI提供两大类型的点对点通信函数：阻塞型(blocking)、非阻塞型(non blocking)。

阻塞型函数需要等待指定操作的实际完成，或至少所涉及的数据已经被MPI系统安全地备份后才返回，如MPI_Send和MPI_Recv，MPI_Send函数返回时表明数据已经发出或者已经被MPI系统复制，随后对发送对象的修改不会影响所发送的数据，而MPI_Recv返回则表明数据已经接收到并且可以立即使用。阻塞型函数的操作是非局部的，它的完成需要与其它进程进行通信。

非阻塞型函数的调用总是立即返回，而实际操作则由MPI系统在后台进行。非阻塞型函数的名称一般为MPI_Ixxx，如MPI_Isend和MPI_Irecv。在调用了一个非阻塞型通信函数后，程序可以去进行其它任务，但必须在之后调用其它函数，如MPI_Wait和MPI_Test等，来等待操作完成或者查询操作的完成情况，在操作完成之前对相关对象的操作是不安全的，因为它们可能与正在后台进行的通信发生冲突。非阻塞型函数的调用是局部的，它不与其它进程发生通信，使用非阻塞型函数可以实现计算与通信重叠。

对于点对点的消息发送，MPI提供4种发送模式，这4种模式的函数具有相同的参数，但它们发送消息的方式不同。每种模式都有阻塞型发送函数和非阻塞性发送函数，如下表所示：

通信模式	阻塞型	非阻塞型
标准模式	MPI_Send	MPI_Isend
缓冲模式	MPI_Bsend	MPI_Ibsend
同步模式	MPI_Ssend	MPI_Issend
就绪模式	MPI_Rsend	MPI_Irsend

而所有模式的阻塞型接收函数和非阻塞型接受函数都各自只有一个：MPI_Recv和MPI_Irecv。不同类型的发送函数的主要区别在于：1、是否使用数据缓冲区；2、函数何时返回。

§阻塞型通信

§1、标准模式(standard mode)

在标准模式下，是否使用数据缓冲区以及对数据缓冲区的管理都是由 MPI 自身决定的，用户无法控制。根据 MPI 是否选择缓存发送数据，可以将发送操作完成的标准可以分为下面两种情况：

MPI 缓存数据 - 在这种情况下，发送操作不管接受操作是否执行，都可以进行，并且发送操作不需要接收操作收到数据就可以成功返回。
MPI 不缓存数据 - 缓存数据是需要付出代价的，它会延长通信的时间，并且缓冲区并不是总能得到的，所以 MPI 可以选择不缓存数据。在这种情况下，只有当接收操作被调用，并且发送的数据完全到达接收缓冲区后，发送操作才算完成。

标准模式对应的阻塞型单点发送函数为：

1
2
3
4
5
6
7
8


int MPI_Send(
    void *buf,              // 发送缓冲区的起始地址
    int count,              // 发送数据的个数
    MPI_Datatype datatype,  // 发送数据的数据类型
    int dest,               // 目标进程号
    int tag,                // 消息标签
    MPI_Comm comm           // 通信器
);

所有模式的阻塞型接收函数都是同一个：

1
2
3
4
5
6
7
8
9


int MPI_Recv(
    void *buf,              // 接收缓冲区
    int count,              // 数据个数上限，具体收到的个数由MPI_Get_count获取
    MPI_Datatype datatype,  // 接收的数据类型
    int source,             // 源进程号
    int tag,                // 消息标签
    MPI_Comm comm,          // 通信器
    MPI_Status *status      // 接收状态，不关心时可以传MPI_STATUS_IGNORE
);

例1、环状发送和接收消息

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34


// 文件名：ring.cpp
#include <mpi.h>
#include <iostream>
#include <random>

int main(int argc, char *argv[]) {
    MPI_Init(&argc, &argv);

    int size = 0, rank = 0;
    MPI_Comm_size(MPI_COMM_WORLD, &size);
    MPI_Comm_rank(MPI_COMM_WORLD, &rank);

    int token = 0;

    if (0 == rank) {
        token = std::random_device{}();
    } else { // 其它进程先接收数据再发送数据
        MPI_Recv(&token, 1, MPI_INT, rank - 1, 0, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
        std::cout << "process " << rank << " received a token " << token <<
                  " from process " << rank - 1 << "\n";
    }

    MPI_Send(&token, 1, MPI_INT, (rank + 1) % size, 0, MPI_COMM_WORLD);

    // 第0个进程接收最后一个进程发送的数据
    if (0 == rank) {
        MPI_Recv(&token, 1, MPI_INT, size - 1, 0, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
        std::cout << "process " << rank << " received a token " << token <<
                  " from process " << size - 1 << "\n";
    }

    MPI_Finalize();
    return 0;
}

编译：mpicxx ring.cpp -o ring

4个进程执行：

1
2
3
4
5


$ mpiexec -n 4 ./ring
process 1 received a token 1502712217 from process 0
process 2 received a token 1502712217 from process 1
process 3 received a token 1502712217 from process 2
process 0 received a token 1502712217 from process 3

使用MPI_Send和MPI_Recv时要特别小心，使用不当很容易引起发送和接收操作不匹配，从而导致程序死锁，比如：

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23


#include <mpi.h>

int main(int argc, char *argv[]) {
    MPI_Init(&argc, &argv);

    int size = 0, rank = 0;
    MPI_Comm_size(MPI_COMM_WORLD, &size);
    MPI_Comm_rank(MPI_COMM_WORLD, &rank);

    constexpr int N = 100;
    int buf1[N];
    int buf2[N];
    
    int dst = (rank + 1) % size;
    int src = (rank - 1 + size) % size;
    int tag = 123;
    
    MPI_Send(buf1, N, MPI_INT, dst, tag, MPI_COMM_WORLD);
    MPI_Recv(buf2, N, MPI_INT, src, tag, MPI_COMM_WORLD, MPI_STATUS_IGNORE);

    MPI_Finalize();
    return 0;
}

如果MPI系统选择缓存数据再发送则不会有问题，如果MPI系统选择直接发送则程序在MPI_Send一行发生死锁，此时所有进程都在等待对方接收数据。如果交换上面程序的MPI_Send和MPI_Recv则在任何情况下都会死锁，此时所有进程都在等待对方发送数据。

§2、缓冲模式(buffered mode)

在缓冲模式下，MPI系统将消息复制到一个用户提供的缓冲区，然后立即返回，消息的发送由MPI系统在后台进行。用户必须确保所提供的缓冲区足以容下发送的消息。缓冲模式发送操作是局部的，因为函数不需要与接收主联络即可返回。如果希望直接对通信缓冲区进行控制，可以使用缓冲通信模式。

缓冲模式的发送函数为MPI_Bsend，参数与MPI_Send完全一致，接收函数仍然为MPI_Recv。在使用MPI_Bsend之前，需要指定缓冲区，否则程序报类似于下面的错误：

Fatal error in MPI_Bsend: Invalid buffer pointer, error stack:
MPI_Bsend(214).......: MPI_Bsend(buf=0x7ffdff7c2d84, count=1, MPI_INT, dest=1, tag=99, MPI_COMM_WORLD) failed
MPIR_Bsend_isend(311): Insufficient space in Bsend buffer; requested 4; total buffer size is 0

指定缓冲区使用MPI_Buffer_attach：

1
2
3
4


int MPI_Buffer_attach(
    void *buffer,   // 缓冲区地址
    int size        // 缓冲区大小（以字节为单位）
);

回收缓冲区使用MPI_Buffer_detach：

1
2
3
4


int MPI_Buffer_detach(
    void **buffer,  // 返回缓冲区地址
    int *size       // 返回缓冲区大小
);

回收操作是阻塞调用，它会一直等到使用该缓存的消息发送完成之后才返回。只有调用返回之后，用户才可以重新使用该缓冲区或者将缓冲区释放。

缓冲区大小的计算比较麻烦，总缓冲区的大小应当至少为所有未完成的MPI_BSend所需缓冲区大小的总和加上通信必须的消息头部大小，头部大小由MPI_BSEND_OVERHEAD给出。所传输的数据大小可以使用MPI_Pack_size获得。

1
2
3
4
5
6


int MPI_Pack_size(
    int incount,            // 数据的个数
    MPI_Datatype datatype,  // 数据的类型
    MPI_Comm comm,          // 通信域
    int *size               // 数据所需要的空间，以字节为单位
);

缓冲区大小计算方式大致为：

1
2
3
4
5
6


MPI_Pack_size(c1, datatype1, comm, &s1);
MPI_Pack_size(c2, datatype2, comm, &s2);
size = s1 + s2 + 2 * MPI_BSEND_OVERHEAD;
// ...
MPI_Bsend( ..., count=c1, datatype1,  ... );
MPI_Bsend( ..., count=c2, datatype2,  ... );

例2：使用MPI_Bsend发送消息

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36


#include <iostream>
#include <random>
#include <mpi.h>

int main(int argc, char *argv[]) {
    MPI_Init(&argc, &argv);

    int rank = 0;
    MPI_Comm_rank(MPI_COMM_WORLD, &rank);
    
    int token = 0;
    
    if(rank == 0) {
        token = std::random_device{}();
        
        // 计算缓冲区大小
        int size;
        MPI_Pack_size(1, MPI_INT, MPI_COMM_WORLD, &size);
        size += MPI_BSEND_OVERHEAD;
        int *tmp_buffer = new int[size];
        
        // 指定缓冲区
        MPI_Buffer_attach(tmp_buffer, size);
        MPI_Bsend(&token, 1, MPI_INT, 1, 123, MPI_COMM_WORLD);
        
        // 回收缓冲区
        MPI_Buffer_detach(&tmp_buffer, &size);
        delete [] tmp_buffer;
    } else if (rank == 1) {
        MPI_Recv(&token, 1, MPI_INT, 0, 123, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
        std::cout << "token received from process 0 is " << token << std::endl;
    }

    MPI_Finalize();
    return 0;
}

运行结果：token received from process 0 is 973762974。

§3、同步模式(synchronous mode)

在标准模式的基础上要求确认接收方已经开始接收数据后函数调用才返回。

§4、就绪模式(ready mode)

调用就绪模式发送消息时必须确保接收方已经处于就绪状态(正在等待接收该消息)，否则会产生一个错误。就绪模式的目的是在一些以同步方式工作的并行系统上，由于发送时可以假设接收方已经准备好接收而减少一些接收开销。就绪通信模式的特殊之处在于它要求接受操作先于发送操作而被启动，因此在一个正确的程序中，一个就绪发送可以被一个标准发送替代，它对程序的正确性没有影响，而对程序的性能有影响。

§非阻塞型通信

§概述

当一个阻塞型通信函数正确返回后，必定有：发送的对象或数组已经被接收或者复制，可以更改该对象或数组；接收的对象或数组已经完整，可以马上使用。与阻塞通信不同，非阻塞通信不必等到通信操作完成完成便可以返回，相对应的通信操作会交给MPI系统在后台进行，发送进程可以去做其它任务，但发送的数组或者对象可能仍然在被后台的发送任务使用，不能对它们进行更改，必须等待发送操作完成。

非阻塞型发送函数也有标准模式、缓冲模式、同步模式和就绪模式四种模式的函数，函数名称分别为MPI_Isend、MPI_Ibsend、MPI_Issend和MPI_Irsend。阻塞型接收函数也只有一个，为MPI_Irecv。点对点通信中，发送类型和接受类型可以不同，阻塞型函数发送的消息可以由非阻塞函数接收，非阻塞型函数发送的消息也可以由阻塞型函数接收。

非阻塞型发送函数的声明为：

1
2
3
4
5
6
7
8
9


MPI_Ixsend(
    void *buf,              // 发送缓冲区起始地址
    int count,              // 发送数据个数
    MPI_Datatype datatype,  // 发送数据的数据类型
    int dest,               // 目标进程号
    int tag,                // 消息标签
    MPI_Comm comm,          // 通信器
    MPI_Request *request   // 返回的非阻塞通信对象
);

非阻塞型接收函数的声明为：

1
2
3
4
5
6
7
8
9


MPI_Irecv(
    void *buf,              // 接受缓冲区的起始地址
    int count,              // 接受数据的最大个数
    MPI_Datatype datatype,  // 数据类型
    int source,             // 源进程号
    int tag,                // 消息标签
    MPI_Comm comm,          // 通信器
    MPI_Request *request    // 非阻塞通信对象
);

§通信完成情况

由于非阻塞通信返回并不意味着该通信已经完成，因此MPI提供了一个非阻塞通信对象MPI_Request来查询通信的状态。通过结合MPI_Request和MPI提供的一些函数，可以检测非阻塞通信的完成情况。

对于单个非阻塞通信来说，可以使用下面两个函数来等待或者检测非阻塞通信。其中MPI_Wait会阻塞当前进程，一直等到相应的非阻塞通信完成之后再返回。MPI_Test用来检测通信是否完成，它会立即返回，不会阻塞当前进程。如果通信完成，将flag置为true，如果通信还没完成，则将flag置为false。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10


MPI_Wait(
    MPI_Request *request,      // 非阻塞通信对象
    MPI_Status *status         // 返回的状态
);

MPI_Test(
    MPI_Request *request,      // 非阻塞通信对象
    int *flag,                 // 操作是否完成，完成 - true，未完成 - false
    MPI_Status *status         // 返回的状态
);

对于多个非阻塞通信，MPI也提供了相应的等待或者检测函数。MPI_Waitany用来等待多个非阻塞通信中的任何一个的完成，MPI_Waitall等待所有的非阻塞通信完成。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12


MPI_Waitany(
    int count,                        // 非阻塞通信对象的个数
    MPI_Request array_of_requests[],  // 非阻塞通信对象数组
    int *index,                       // 通信完成的对象在数组中的索引
    MPI_Status status[]               // 返回的状态
);

MPI_Waitall(
    int count，                       // 非阻塞通信对象的个数
    MPI_Request array_of_requests[],  // 非阻塞通信对象数组
    MPI_Status array_of_status[]      // 状态数组
);

相应的通信检测函数MPI_Testany、MPI_Testall、MPI_Probe以及通信取消函数MPI_Cancel等可以查看文档。

§例3：自适应区间积分

之前我们使用定积分的定义来近似计算$\pi$的值：

$$ \int_0^1\cfrac{4}{1 + x^2}dx = \lim_{n \rightarrow +\infty} h\sum_{i=0}^{n-1} f(\delta_i) $$ 其中$h=1/n$为积分步长；$\delta_i=(i+0.5)h$为区间中点；$f(x)=\cfrac{4}{1 + x^2}$为被积函数。这里将区间[0, 1]等分为n份，每份的上的定积分近似看作一个长方形的面积，长方形的底为$h$，高为$f(\delta_i)$。

此外，对每个小区间$[x_i,x_{i+1}]$，还可以将定积分近似看作一个梯形的面积，此梯形的面积为$(x_{i+1}-x_i)\cfrac{f(x_i)+f(x_{i+1})}{2}$，将所有小区间上的面积相加也可以得到定积分的近似值： $$ \int_a^b f(x) dx = \frac{h}{2} \left( f(a) + f(b) + 2\sum_{i=1}^{n-1} f(x_i) \right) $$ 其中$h=(b-a)/n$为积分步长，$x_i=a+ih$为第i个小区间左端点。

这类方法的问题在于某些小区间划分得太细，而某些小区间的误差又太大，因此有必须采用一种自适应的划分方法。二分法计算积分$\int_a^b f(x) dx$描述如下：

(0) 当前计算区间为$(x_0,x_1)=(a,b)$。

(1) 用梯形面积公式$(x_{i+1}-x_i)\cfrac{f(x_i)+f(x_{i+1})}{2}$计算当前区间上的积分近似值。

(2) 判断近似值的误差，如果误差满足要求，则该值为该区间上的积分近似值，结束当前区间上的计算。

(3) 否则将区间一分为二，令$x_c=(x_0+x_1)/2$，分别对子区间$(x_0,x_c)$和$(x_c,x_1)$重复步骤(1)-(3)，然后把它们的积分近似结果相加作为区间$(x_0,x_1)$上的计算结果。

第(2)步中的误差判断方法为：$v_0-v<3h\varepsilon/(b-a)$。其中$v_0=\cfrac{h}{2}(f(x_0)+f(x_1))$，$v=(v_0+hf(x_c))/2$是下一步的积分近似值，$\varepsilon$为整个积分的误差要求，$h=x_1-x_0$为当前区间积分步长。(参考http://math.ecnu.edu.cn/~jypan/Teaching/SciComp/lect01B.pdf和https://en.wikipedia.org/wiki/Trapezoidal_rule)

串行代码如下：

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70


#include <iostream>
#include <functional>
#include <cmath>

struct Interval {
    double a, b, fa, fb, fc;
    Interval *last, *next;
};

double Integrate(const std::function<double(double)> &f,
                 double a, double b,
                 double eps = 1e-15) {
    eps *= 3 * (b - a);     // 精度计算中的公因式预先算上

    double result = 0;
    Interval *root = new Interval{a, b, f(a), f(b), 0, nullptr, nullptr};
    while (root != nullptr) {
        // 先计算各区间中点函数值
        for (Interval *i = root; i != nullptr; i = i->next)
            i->fc = f(0.5 * (i->a + i->b));

        // 计算各区间积分值，如果精度不够就继续细分
        Interval *i = root;
        while (i != nullptr) {
            double h = i->b - i->a;
            double xc = 0.5 * (i->a + i->b);
            double v0 = 0.5 * h * (i->fa + i->fb);
            double v1 = 0.5 * (v0 + i->fc * h);

            // 误差满足要求或者区间长度太小则计算此区间内的积分并将此区间移出C
            if (fabs(v1 - v0) < h * eps || xc == i->a || xc == i->b) {
                result += v1;
                if (i->next != nullptr)
                    i->next->last = i->last;
                if (i->last != nullptr)
                    i->last->next = i->next;
                Interval *temp = i;
                i = i->next;
                if (temp == root)
                    root = i;
                delete temp;
                continue;
            }

            // 将i拆分成两个区间，i保留左半区间，i后面插入右半区间
            Interval *j = new Interval;
            j->a = xc;
            j->b = i->b;
            j->fa = i->fc;
            j->fb = i->fb;
            i->b = xc;
            i->fb = i->fc;
            j->next = i->next;
            if (i->next != nullptr)
                i->next->last = j;
            i->next = j;
            j->last = i;
            i = j->next;    // 跳过刚插入的j
        }
    }

    return result;
}

int main() {
    double pi = Integrate([](double x) { return 4 / (1 + x*x); }, 0, 1, 1e-15);
    std::cout.precision(16);
    std::cout << "pi is " << pi << "\n";
    return 0;
}

MPI代码如下：

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145


#include <iostream>
#include <functional>
#include <cmath>
#include <cstdlib>
#include <mpi.h>

struct Interval {
    double a, b, fa, fb, fc;
    Interval *last, *next;
};

double Integrate(const std::function<double(double)> &f,
                 double a, double b,
                 double eps = 1e-15) {
    // 散发求值任务，使多次求函数值的过程并行进行。F的功能可以看作请求求值，结果放入v中，当
    // v为nullptr时，功能为等待所有已提交的求值请求完成。
    int np = 0;
    MPI_Comm_size(MPI_COMM_WORLD, &np);
    MPI_Request *request = nullptr;
    MPI_Status *status = nullptr;
    int size = 0;                // 已分配的单元个数
    std::function<void(double, double *)> F =
            [&np, &request, &status, &size](double x, double *v) -> void {
        static int n = 0;        // 目前有多少个未计算函数值的x
        static int count = 0;    // 全局计算请求总数
        static int slave;        // 根据计算请求的总数计算要将计算任务发到哪个进程上

        if (v == nullptr && n > 0) {
            MPI_Waitall(n + n, request, status);
            n = 0;
            return;
        }

        // 需要增加内存
        if (n >= size) {
            size += 128;        // 增加128个单元的内存
            request = (MPI_Request *)realloc(request, 2 * size * sizeof(MPI_Request));
            status = (MPI_Status *)realloc(status, 2 * size * sizeof(MPI_Status));
        }

        slave = count % (np - 1) + 1;   // 不发到0号进程上
        MPI_Isend(&x, 1, MPI_DOUBLE, slave, 1, MPI_COMM_WORLD, request + 2 * n);
        MPI_Irecv(v, 1, MPI_DOUBLE, slave, 1, MPI_COMM_WORLD, request + 2 * n + 1);
        ++n;
        ++count;
    };

    eps *= 3 * (b - a);     // 精度计算中的公因式预先算上
    double result = 0;
    Interval *root = new Interval{a, b, f(a), f(b), 0, nullptr, nullptr};
    while (root != nullptr) {
        // 先计算各区间中点函数值
        for (Interval *i = root; i != nullptr; i = i->next)
            F(0.5 * (i->a + i->b), &i->fc);
        F(0, nullptr);    // 等待所有函数值计算完成

        // 计算各区间积分值，如果精度不够就继续细分
        Interval *i = root;
        while (i != nullptr) {
            double h = i->b - i->a;
            double xc = 0.5 * (i->a + i->b);
            double v0 = 0.5 * h * (i->fa + i->fb);
            double v1 = 0.5 * (v0 + i->fc * h);

            // 误差满足要求或者区间长度太小则计算此区间内的积分并将此区间移出C
            if (fabs(v1 - v0) < h * eps || xc == i->a || xc == i->b) {
                result += v1;
                if (i->next != nullptr)
                    i->next->last = i->last;
                if (i->last != nullptr)
                    i->last->next = i->next;
                Interval *temp = i;
                i = i->next;
                if (temp == root)
                    root = i;
                delete temp;
                continue;
            }

            // 将i拆分成两个区间，i保留左半区间，i后面插入右半区间
            Interval *j = new Interval;
            j->a = xc;
            j->b = i->b;
            j->fa = i->fc;
            j->fb = i->fb;
            i->b = xc;
            i->fb = i->fc;
            j->next = i->next;
            if (i->next != nullptr)
                i->next->last = j;
            i->next = j;
            j->last = i;
            i = j->next;    // 跳过刚插入的j
        }
    }

    delete[] request;
    delete[] status;
    return result;
}

int main(int argc, char *argv[]) {
    auto f = [](double x) { return 4 / (1 + x * x); };  // 被积函数
    double a = 0, b = 1;                                // 积分区间
    double eps = 1e-7;                                  // 积分精度

    int np = 0, id = 0;
    MPI_Init(&argc, &argv);
    MPI_Comm_size(MPI_COMM_WORLD, &np);
    MPI_Comm_rank(MPI_COMM_WORLD, &id);

    // 主从模式至少需要两个进程
    if (np < 2) {
        std::cerr << "This program needs at least two processes!" << std::endl;
        MPI_Abort(MPI_COMM_WORLD, -1);
    }

    // 0号进程作为主进程
    if (0 == id) {
        double pi = Integrate(f, a, b, eps);
        std::cout.precision(16);
        std::cout << "pi is " << pi << ".\n";

        // 通知其它进程退出
        for (int i = 1; i < np; ++i)
            MPI_Send(&pi, 1, MPI_DOUBLE, i, 0, MPI_COMM_WORLD);
    } else {
        while (true) {
            double d = 0;
            MPI_Status status;
            MPI_Recv(&d, 1, MPI_DOUBLE, 0, MPI_ANY_TAG, MPI_COMM_WORLD, &status);
            
            // tag为0时是退出信号，tag为1时是计算任务
            if (status.MPI_TAG == 0)
                break;
            
            // 计算任务时将指定点的函数值求出，然后发送给0号进程
            d = f(d);
            MPI_Send(&d, 1, MPI_DOUBLE, 0, 1, MPI_COMM_WORLD);
        }
    }

    MPI_Finalize();
    return 0;
}

这是典型的主从模式并行计算程序，0号进行负责计算误差、发布计算任务等工作，其它进程为从进程，负责计算各点的函数值。

§文档

[1] https://www.open-mpi.org/doc/current/.

[2] https://www.mpich.org/static/docs/latest/.

[3] https://www.mpi-forum.org/docs/mpi-3.1/mpi31-report.pdf.

MPI教程2

相关文章：