多线程并行快速求解Pi的方法

摘 要 本文首先分别应用不同的求Pi方法分析讲解了WinAPI、OpenMP、MPI三大并行算法中常遇到的一些问题，根据问题提出了相应的解决方案。另外实现了对BBP算法的初步并行化，大大缩减了该算法的运行时间。文章最后利用三种方法求解了蒙特卡洛算法求Pi，并对比分析了三种算法的优缺点。

【关键词】WinAPI OpenMP MPI

1 问题背景与提出

随着技术的发展，单片机越来越难满足人们对于大量数据处理的需求，因此，人们越来越依赖于利用并行计算技术来解决程序规模庞大，运算时间长及数据量大的课题。本文即对当下比较常用的几种并行技术WinAPI、OpenMP、MPI三种并行模式进行讨论课研究，并以计算π值为例，将并行模式与串行模式进行对比，研究并行计算的机理、优缺点及一些常见问题。

2 模型建立

2.1 蒙特卡_思想，蒲风投针实验

（1）取白纸一张，在上面画许多间距为d的等距平行线；

（2）取一根长为l（l

（3）直线与针相交概率p的近似值可用m/n得到，进而可得到圆周率的近似值为2nl/md。

2.2 级数方法

2.3 并行BBP算法

当今世界在进行计算机性能测试时往往选择在一定时间内计算Pi值得位数，而最常用的算法之一就是BBP算法。该算法不需要多精度浮点算术运算的支持，在支持IEEE浮点运算的通用计算机上即可进行计算而且该算法只需要非常少的内存。BBP算法也是目前算法中非常适用于并行计算的算法。

公式为：

3 算法描述与实现

3.1 API实现

Windows实现多线程并行计算时会产生数据竞争，数据竞争（Data racing）导致计算结果不准确。产生错误的原因在于两个线程在同时访问同一内存区域时，且一个线程在进行写操作。为此采用同步技术，利用临界区解决此问题。由于本次试验中，数据计算量较小，所以基本不会产生数据竞争的错误，但是一旦运行计算量较大或者在共享量上运算时间较长的程序时，就会产生数据竞争，在本例中，我们以cout输出程序为引子，观察使用临界区和不使用临界区的区别：

代码如下：

}

pi= count*4.0 / numSteps；

//EnterCriticalSection（&cs）；

cout

piSum+= pi；

// LeaveCriticalSection（&cs）；

我们将输出程序置于piSum值之前，由于cout运行时，在屏幕上显示需要消耗较大时间，所以就会有机会产生数据重叠或者两者读到了一样的piSum值。（如图1所示）。

所以，当我们使用临界区后结果为：（如图2所示）。

另外，临界区的设定也直接影响计算速度，所以尽量缩小临界区有利于提高运行速度。

线程取值过大的影响：

3.2 OpenMP实现级数方法

OpenMP是一种面向共享内存以及分布共享内存的多处理器多线程并行编程语言，是一种能够被应用于显式指导多线程、共享内存并行的应用程序编程接口（如图3所示）。OpenMP具有良好可移植性，支持多种编程语言。parallel for 循环并行化利用编译指导语句parallel for 并行原理是采用工作分配的执行方式，将循环所需要的工作量按一定方式分配到各个线程，所有线程执行工作总和是原串行完成的工作量。parallel for 根据CUP的大小（线程数多少）按工作量平均分配，对于线程数为8个，这里蒙特卡洛方法工作量取numSteps=10，000，000步，则每个线程要计算（1/8）*numSteps步的工作量。要特别注意的是在这里必须避免数据竞争（Data racing），采用的方法是对竞争变量私有化，使用private（t，fun，...）。

核心代码：

double Pi（double x）

{

omp_set_num_threads（numThreads）；

int sign = 1；

#pragma omp parallel for reduction（+：c[omp_get_thread_num（）]）

for （int i = 2； i < N； i++）

{

c[omp_get_thread_num（）] += 4 * x；

sign = sign * （-1）；

x = sign /double （2 * i - 1）；

}

这个算法中，我们发现运行后，串行时间和并行时间是基本一致的：

这个问题在于x，sign的值为共享量，当不同的线程在调用时存在等待的时间，所以如果将这些值私有化，则可以解决这些问题。由于使用private时，各个私有变量必须在for中有定义，即在进入并行区时，for循环中的所有私有变量是需要重新定义的，因此我们必须使得x，sign等变量在定义时全部使用已知量来定义

#pragma omp parallel for private（x，y）reduction（+：pi）

for （int i = 1； i < N； i += 2）

{

y = 1 / double（2 * i - 1）；