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

于朋

摘 要 本文首先分別應(yīng)用不同的求Pi方法分析講解了WinAPI、OpenMP、MPI三大并行算法中常遇到的一些問題，根據(jù)問題提出了相應(yīng)的解決方案。另外實(shí)現(xiàn)了對BBP算法的初步并行化，大大縮減了該算法的運(yùn)行時間。文章最后利用三種方法求解了蒙特卡洛算法求Pi，并對比分析了三種算法的優(yōu)缺點(diǎn)。

【關(guān)鍵詞】WinAPI OpenMP MPI

1 問題背景與提出

隨著技術(shù)的發(fā)展，單片機(jī)越來越難滿足人們對于大量數(shù)據(jù)處理的需求，因此，人們越來越依賴于利用并行計算技術(shù)來解決程序規(guī)模龐大，運(yùn)算時間長及數(shù)據(jù)量大的課題。本文即對當(dāng)下比較常用的幾種并行技術(shù)WinAPI、OpenMP、MPI三種并行模式進(jìn)行討論課研究，并以計算π值為例，將并行模式與串行模式進(jìn)行對比，研究并行計算的機(jī)理、優(yōu)缺點(diǎn)及一些常見問題。

2 模型建立

2.1 蒙特卡羅思想，蒲風(fēng)投針實(shí)驗

（1）取白紙一張，在上面畫許多間距為d的等距平行線；

（2）取一根長為l（l

（3）直線與針相交概率p的近似值可用m/n得到，進(jìn)而可得到圓周率的近似值為2nl/md。

2.2 級數(shù)方法

2.3 并行BBP算法

當(dāng)今世界在進(jìn)行計算機(jī)性能測試時往往選擇在一定時間內(nèi)計算Pi值得位數(shù)，而最常用的算法之一就是BBP算法。該算法不需要多精度浮點(diǎn)算術(shù)運(yùn)算的支持，在支持IEEE浮點(diǎn)運(yùn)算的通用計算機(jī)上即可進(jìn)行計算而且該算法只需要非常少的內(nèi)存。BBP算法也是目前算法中非常適用于并行計算的算法。

公式為：

3 算法描述與實(shí)現(xiàn)

3.1 API實(shí)現(xiàn)

Windows實(shí)現(xiàn)多線程并行計算時會產(chǎn)生數(shù)據(jù)競爭，數(shù)據(jù)競爭（Data racing）導(dǎo)致計算結(jié)果不準(zhǔn)確。產(chǎn)生錯誤的原因在于兩個線程在同時訪問同一內(nèi)存區(qū)域時，且一個線程在進(jìn)行寫操作。為此采用同步技術(shù)，利用臨界區(qū)解決此問題。由于本次試驗中，數(shù)據(jù)計算量較小，所以基本不會產(chǎn)生數(shù)據(jù)競爭的錯誤，但是一旦運(yùn)行計算量較大或者在共享量上運(yùn)算時間較長的程序時，就會產(chǎn)生數(shù)據(jù)競爭，在本例中，我們以cout輸出程序為引子，觀察使用臨界區(qū)和不使用臨界區(qū)的區(qū)別：

代碼如下：

}

pi= count*4.0 / numSteps；

//EnterCriticalSection（&cs）；

cout << piSum<< " " << threadID << endl；

piSum+= pi；

// LeaveCriticalSection（&cs）；

我們將輸出程序置于piSum值之前，由于cout運(yùn)行時，在屏幕上顯示需要消耗較大時間，所以就會有機(jī)會產(chǎn)生數(shù)據(jù)重疊或者兩者讀到了一樣的piSum值。（如圖1所示）。

所以，當(dāng)我們使用臨界區(qū)后結(jié)果為：（如圖2所示）。

另外，臨界區(qū)的設(shè)定也直接影響計算速度，所以盡量縮小臨界區(qū)有利于提高運(yùn)行速度。

線程取值過大的影響：

3.2 OpenMP實(shí)現(xiàn)級數(shù)方法

OpenMP是一種面向共享內(nèi)存以及分布共享內(nèi)存的多處理器多線程并行編程語言，是一種能夠被應(yīng)用于顯式指導(dǎo)多線程、共享內(nèi)存并行的應(yīng)用程序編程接口（如圖3所示）。OpenMP具有良好可移植性，支持多種編程語言。parallel for 循環(huán)并行化利用編譯指導(dǎo)語句parallel for 并行原理是采用工作分配的執(zhí)行方式，將循環(huán)所需要的工作量按一定方式分配到各個線程，所有線程執(zhí)行工作總和是原串行完成的工作量。parallel for 根據(jù)CUP的大小（線程數(shù)多少）按工作量平均分配，對于線程數(shù)為8個，這里蒙特卡洛方法工作量取numSteps=10，000，000步，則每個線程要計算（1/8）*numSteps步的工作量。要特別注意的是在這里必須避免數(shù)據(jù)競爭（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）；

}

這個算法中，我們發(fā)現(xiàn)運(yùn)行后，串行時間和并行時間是基本一致的：

這個問題在于x，sign的值為共享量，當(dāng)不同的線程在調(diào)用時存在等待的時間，所以如果將這些值私有化，則可以解決這些問題。由于使用private時，各個私有變量必須在for中有定義，即在進(jìn)入并行區(qū)時，for循環(huán)中的所有私有變量是需要重新定義的，因此我們必須使得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）；

pi = pi + 4 * y；

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

pi = pi + 4 * x；

}

因此程序完美的解決了并行問題。（如表1所示）。

3.3 BBP算法的實(shí)現(xiàn)

MPI是專門為集群和多節(jié)點(diǎn)并行計算語言，運(yùn)行效率高，實(shí)現(xiàn)方式多樣，可以進(jìn)行主從式、并列式以及流水線式等方式的實(shí)現(xiàn)。在利用計算機(jī)多核CPU，模擬多個節(jié)點(diǎn)的實(shí)現(xiàn)方式。改造成主從式程序，利用0號節(jié)點(diǎn)作為主節(jié)點(diǎn)收發(fā)數(shù)據(jù)，也參與計算，而其他節(jié)點(diǎn)只進(jìn)行計算，為負(fù)載均衡選擇詳盡的計算量運(yùn)行到不同的節(jié)點(diǎn)上。

根據(jù)BBP算法（如圖4所示）我們將整個運(yùn)算分成4個部分，分別用一個線程進(jìn)行運(yùn)算，其運(yùn)算總時間約等于線程中運(yùn)行時間最長的那個（計算結(jié)果如表2所示）。

3.4 API、OpenMP、MPI對比實(shí)驗

我們以蒙塔卡羅算法為例，分別采用API、OpenMP、MPI三種并行方式對值進(jìn)行計算，來對比其運(yùn)算效率和運(yùn)算精度（計算結(jié)果如表3所示）。

4 結(jié)論

（1）MPI模擬多節(jié)點(diǎn)計算的速度是最快的，而且計算結(jié)果也是最接近穿行結(jié)果的，所以MPI適合計算大規(guī)模的較為精確的求解問題。

（2）WinAPI實(shí)現(xiàn)用臨界區(qū)的效果最差，而且WinAPI的計算速度很大程度也和臨界區(qū)的設(shè)定有關(guān)，盡量縮小臨界區(qū)有利于提高并行速度。

（3）WinAPI實(shí)現(xiàn)用線程號為每個線程分配不同的任務(wù)，分配過程需要人為干預(yù)，對于函數(shù)復(fù)雜的程序來說，實(shí)現(xiàn)繁瑣，不利于使用。

（4）數(shù)值計算時要根據(jù)實(shí)際情況選擇和改造算法。比如在計算值就比e更適合并行。而且每種并行方法都有它的特使要求，比如在使用parallel for private時，由于變量進(jìn)入for循環(huán)后屬于重新定義，所以不能出現(xiàn)自己為自己賦值的情況，需要按照程序一步步展開來寫，相對繁瑣。

參考文獻(xiàn)

[1]多核系列教材編寫組.多核程序設(shè)計[M].北京：清華大學(xué)出版社，2007.

[2]朱建偉，劉榮.多線程并行快速求解e值的六種方法[J].現(xiàn)代計算機(jī)，2013.

[3]張翔圓.周率Pi的BBP多核并行算法實(shí)現(xiàn)[J].普洱學(xué)院學(xué)報，2013.

作者單位

中國石油大學(xué)（華東）理學(xué)院 山東省青島市 266500