题目描述
该程序使用 MPI (Message Passing Interface) 并行计算 $\pi$ 的近似值。它基于一个经典的微积分公式：
$$ \pi = \int_0^1 \frac{4}{1+x^2} dx $$
程序通过数值积分的方法（具体来说是中点矩形法）来逼近这个定积分的值。它将积分区间 $[0, 1]$ 分成 $n$ 个小区间，然后将这些小区间的计算任务分配给多个 MPI 进程并行处理，最后汇总结果得到 $\pi$ 的近似值。
算法与实现解析
(基于MPI的并行中点矩形法数值积分)
该算法的核心思想是将计算密集型的求和任务分解到多个 MPI 进程上，以加速 $\pi$ 值的计算。
1. 数学原理：定积分求 $\pi$
我们知道反正切函数 $\arctan(x)$ 的导数是 $\frac{1}{1+x^2}$。因此：
$$ \int \frac{1}{1+x^2} dx = \arctan(x) + C $$
所以，定积分：
$$ \int_0^1 \frac{4}{1+x^2} dx = 4 \left[ \arctan(x) \right]_0^1 = 4 (\arctan(1) - \arctan(0)) = 4 \left( \frac{\pi}{4} - 0 \right) = \pi $$
这就是程序计算 $\pi$ 的数学基础。
2. 数值积分方法：中点矩形法
为了在计算机上计算这个积分，我们使用数值积分。中点矩形法是一种常用的方法。
其基本思想是将积分区间 $[a, b]$ 分成 $n$ 个等宽的小区间，每个区间的宽度为 $h = (b-a)/n$。对于第 $i$ 个小区间 $[x_{i-1}, x_i]$，我们取该区间的中点 $x_i^* = (x_{i-1} + x_i)/2$ 处函数值 $f(x_i^*)$ 与区间宽度 $h$ 的乘积，作为该小区间上积分的近似值。
总的积分近似值为所有这些小矩形面积之和：
$$ \int_a^b f(x) dx \approx \sum_{i=1}^{n} f(x_i^*) h $$
在本题中，$a=0, b=1, f(x) = \frac{4}{1+x^2}$。
所以，$h = (1-0)/n = 1/n$。
第 $i$ 个子区间为 $[(i-1)h, ih]$ (其中 $i$ 从 $1$ 到 $n$)。
该子区间的中点为 $x_i^* = (i-1)h + \frac{h}{2} = ih - \frac{h}{2} = h(i - 0.5)$。
因此，$\pi$ 的近似计算公式为：
$$ \pi \approx \sum_{i=1}^{n} \frac{4}{1 + (h(i-0.5))^2} \cdot h $$
$$ \pi \approx h \sum_{i=1}^{n} \frac{4}{1 + (h(i-0.5))^2} $$
3. MPI 并行化思路
MPI (Message Passing Interface) 是一种标准化的并行编程库，允许在分布式内存系统中的多个进程之间进行通信。
-
初始化 (Initialization):
MPI_Init(&argc, &argv);: 初始化 MPI 执行环境。MPI_Comm_rank(comm, &rank);: 获取当前进程在通信子comm(这里是MPI_COMM_WORLD，包含所有进程的默认通信子) 中的秩 (rank)，即进程的唯一标识符，从 0 到size-1。MPI_Comm_size(comm, &size);: 获取通信子comm中的总进程数。
-
参数广播 (Broadcast):
if (rank == root) { n = 2e9; }: 根进程 (rank 0) 设置积分区间的划分数量 $n$。注意这里 $n$ 是一个int，而2e9(即 $2 \times 10^9$) 是一个很大的整数。MPI_Bcast(&n, 1, MPI_INT, root, comm);: 根进程 (rankroot) 将变量n的值广播给通信子comm中的所有其他进程。这样，每个进程都知道了总的划分数量 $n$。
-
任务分配与局部计算 (Work Distribution & Local Computation):
double h = 1.0 / n;: 每个进程根据接收到的 $n$ 计算子区间的宽度 $h$。double sum = 0.0;: 每个进程初始化自己的局部和sum。-
The core loop:
c++ for (int i = rank + 1; i <= n; i += size) { double x = h * (double(i) - 0.5); sum += 4.0 / (1 + x * x); }
这是并行计算的关键。循环变量i代表第 $i$ 个子区间 (从 1 到 $n$)。
每个进程rank只计算一部分子区间的和。具体来说：- 进程 0 计算 $i = 1, 1+size, 1+2size, \dots$
- 进程 1 计算 $i = 2, 2+size, 2+2size, \dots$
- …
- 进程
rank计算 $i = rank+1, rank+1+size, rank+1+2size, \dots$ - …
- 进程
size-1计算 $i = size, 2size, 3size, \dots$ (因为rank+1对于size-1就是size)
这种分配方式称为循环或交错分配 (cyclic/interleaved distribution)。它确保了如果 $n$ 不是
size的整数倍时，任务也能相对均匀地分配。
double x = h * (double(i) - 0.5);: 计算第 $i$ 个子区间的中点 $x_i^*$。
sum += 4.0 / (1 + x * x);: 累加当前进程负责的 $f(x_i^*)$ 的值。 -
double mypi = sum * h;: 每个进程将自己计算的局部和sum乘以 $h$，得到该进程贡献的 $\pi$ 的一部分，记为mypi。
-
结果汇总 (Reduction):
MPI_Reduce(&mypi, &pi, 1, MPI_DOUBLE, MPI_SUM, root, comm);:
这是一个归约操作。&mypi: 每个进程提供的输入数据的地址（即各自的局部 $\pi$ 值）。&pi: 根进程root(rank 0) 接收归约结果的地址。其他进程此参数无意义，但仍需提供。1: 要归约的数据元素个数（这里是一个double值）。MPI_DOUBLE: 归约数据的 MPI 数据类型。MPI_SUM: 归约操作的类型，这里是求和。所有进程的mypi值将被加起来。root: 接收最终结果的进程的秩。comm: 执行归约操作的通信子。
执行后，根进程的变量pi中将存储所有进程mypi的总和，即最终的 $\pi$ 近似值。
-
输出与终止 (Output & Finalization):
if (rank == 0) { printf("pi is approximately %.16f\n", pi); }: 根进程 (rank 0) 打印计算得到的 $\pi$ 值。MPI_Finalize();: 结束 MPI 执行环境，进行必要的清理工作。
时间复杂度
假设有 size 个进程，总共要计算 $n$ 个点。
每个进程大致处理 $n/size$ 个点。每个点的计算是常数时间 $O(1)$。
所以，并行计算部分的时间复杂度是 $O(n/size)$。
MPI_Bcast 和 MPI_Reduce 的时间复杂度通常与进程数 size 有关，对于一个较小的数据量，可以认为是 $O(\log size)$ 或 $O(size)$，具体取决于 MPI 实现和底层网络。
因此，总的时间复杂度主要由计算部分决定，即 $O(n/size)$，体现了并行计算的加速效果。
C++ 代码
#include <mpi.h>
#include <stdio.h>
#include <stdlib.h>
int main(int argc, char* argv[]) {
int rank, size;
int data_to_bcast; // 这个变量在您的代码中声明了但未使用，n 被用来广播
int root = 0;
MPI_Comm comm = MPI_COMM_WORLD;
int n; // 积分区间的划分数量
MPI_Init(&argc, &argv); // 初始化MPI环境
MPI_Comm_rank(comm, &rank); // 获取当前进程的秩
MPI_Comm_size(comm, &size); // 获取总进程数
// 根进程 (rank 0) 初始化 n
if (rank == root) {
n = 2e9; // 设定一个较大的 n 以获得较精确的结果
}
// 根进程将 n 的值广播给所有其他进程
// 参数: 数据地址, 数据个数, 数据类型, 根进程秩, 通信子
MPI_Bcast(&n, 1, MPI_INT, root, comm);
// 每个进程计算子区间的宽度 h
double h = 1.0 / n;
double sum = 0.0; // 每个进程的局部和
// 循环计算各自负责的区间的函数值之和
// i 从 rank+1 开始，每次增加 size，确保任务分配
// 例如，size=4:
// rank 0: i = 1, 5, 9, ...
// rank 1: i = 2, 6, 10, ...
// rank 2: i = 3, 7, 11, ...
// rank 3: i = 4, 8, 12, ...
for (int i = rank + 1; i <= n; i += size) {
double x = h * (double(i) - 0.5); // 计算第i个子区间的中点
sum += 4.0 / (1 + x * x); // 累加 f(x_midpoint)
}
// 每个进程计算自己那部分的 pi 值 (sum * h)
double mypi = sum * h;
double pi; // 用于存放最终的 pi 结果 (仅在根进程有意义)
// 将所有进程的 mypi 值相加，结果存储在根进程的 pi 变量中
// 参数: 发送缓冲地址, 接收缓冲地址, 数据个数, 数据类型, 操作类型, 根进程秩, 通信子
MPI_Reduce(&mypi, &pi, 1, MPI_DOUBLE, MPI_SUM, root, comm);
// 根进程打印结果
if (rank == 0) {
printf("pi is approximately %.16f\n", pi);
}
MPI_Finalize(); // 结束MPI环境
return 0;
}
