蔡曉龍（洛陽師范學(xué)院，河南　洛陽　471022）

基于MIC的RAN2并行化設(shè)計(jì)與實(shí)現(xiàn)

蔡曉龍
（洛陽師范學(xué)院，河南洛陽471022）

摘要：RAN2是線性同余隨機(jī)數(shù)發(fā)生器中的一種，因其周期長且隨機(jī)性好廣為人們應(yīng)用.為了提升RAN2的產(chǎn)生速度，本文在研究RAN2串行算法的基礎(chǔ)上，利用Simple skip ahead方法對其進(jìn)行并行化.實(shí)驗(yàn)結(jié)果顯示，并行化后的RAN2通過了TestU01的455項(xiàng)測試，與串行結(jié)果相同.相對于CPU單線程，MIC平臺下的最高加速比為17.25.

關(guān)鍵詞：隨機(jī)數(shù)發(fā)生器；RAN2；并行化；MIC

能夠產(chǎn)生隨機(jī)數(shù)序列的軟件或硬件或者二者的結(jié)合被稱為隨機(jī)數(shù)發(fā)生器RNG（random number generator)[1，2].LCG是目前主流隨機(jī)數(shù)發(fā)生器中的一種，其優(yōu)點(diǎn)是隨機(jī)性好、產(chǎn)生速度快等，但周期較短[3].因此，如何提高隨機(jī)數(shù)發(fā)生器的品質(zhì)及其發(fā)生速度成了人們首要關(guān)心的問題，經(jīng)過不斷研究，人們提出了組合隨機(jī)數(shù)發(fā)生器的思想[4，5].Press和Teukolsky在1992年提出了RAN2，這種隨機(jī)數(shù)發(fā)生器利用兩個(gè)隨機(jī)數(shù)序列相加，并用其中一個(gè)隨機(jī)數(shù)發(fā)生器的模作為模[6].然而由于RAN2算法復(fù)雜度的增加，使得其產(chǎn)生隨機(jī)數(shù)的速率較慢.

現(xiàn)今計(jì)算機(jī)的發(fā)展可以明顯分為兩個(gè)階段:串行計(jì)算時(shí)代和并行計(jì)算時(shí)代[7].隨著基于多核處理器結(jié)構(gòu)的計(jì)算機(jī)成為市場的主流，并行計(jì)算將會是計(jì)算機(jī)科學(xué)未來發(fā)展極其重要的一個(gè)方向，在提高計(jì)算機(jī)硬件資源利用率的同時(shí)，可顯著的改善算法的性能、程序的執(zhí)行速度等等.近年來，并行化技術(shù)已經(jīng)漸漸開始應(yīng)用于隨機(jī)數(shù)發(fā)生器領(lǐng)域[8，9]，Bradley T，du Toit J和Giles M等人在2010年將隨機(jī)數(shù)發(fā)生器的并行化方法歸納為Simple skip ahead，Strided skip ahead 和Hybrid三種[10].但針對每一類隨機(jī)數(shù)發(fā)生器具體的使用方法、并行化實(shí)施方案、并行化效果并沒有進(jìn)行闡述和說明.

本文就RAN2運(yùn)行速度緩慢的問題，在前人工作的基礎(chǔ)上利用Simple skip ahead方法實(shí)現(xiàn)了RAN2的并行化，并首次基于MIC進(jìn)行了性能測試分析，主要思想是通過在一個(gè)周期內(nèi)劃分線程和分配任務(wù)，每個(gè)線程獨(dú)立的產(chǎn)生一段周期內(nèi)的隨機(jī)數(shù)的子序列達(dá)到并行化.本文研究的主要內(nèi)容包括對于RAN2串行算法的研究和并行算法的設(shè)計(jì)；CPU平臺下并行程序的開發(fā)及性能測試；MIC平臺上的移植及性能測試；通過目前最為完備的隨機(jī)數(shù)發(fā)生器測試庫TestU01進(jìn)行隨機(jī)數(shù)性能測試，以此來保證并檢驗(yàn)并行化的過程的正確性；最終實(shí)驗(yàn)數(shù)據(jù)表明，相對于串行的RAN2隨機(jī)數(shù)發(fā)生器，并行后的RAN2隨機(jī)數(shù)發(fā)生器速度有了明顯的提升.

1 RAN2串行算法及并行化方法

1.1 RAN2串行算法

RAN2隨機(jī)數(shù)發(fā)生器是由2個(gè)相互獨(dú)立的LCG隨機(jī)數(shù)發(fā)生器組合而成的，通過將2個(gè)LCG發(fā)生器的結(jié)果進(jìn)行合并就可以得到最終的RAN2隨機(jī)數(shù).RAN2隨機(jī)數(shù)發(fā)生器的遞推公式如下:

參數(shù)為:

a1=40014，m1=2147483563; a2=40692，m2=2147483399;

RAN2隨機(jī)數(shù)發(fā)生器的周期為兩個(gè)隨機(jī)數(shù)發(fā)生器周期的最小公倍數(shù)：T=2.3×1018，通過給定初始化種子X0，Y0就可以利用以上遞推公式不斷生成隨機(jī)數(shù).當(dāng)產(chǎn)生大量隨機(jī)數(shù)時(shí)，由于初始化時(shí)間很短故可忽略不計(jì)，假設(shè)隨機(jī)數(shù)發(fā)生器運(yùn)行一次的時(shí)間為單位1，則理論上產(chǎn)生random_num個(gè)隨機(jī)數(shù)所需時(shí)間約為O（random_num).

1.2 RAN2并行化實(shí)現(xiàn)

圖1 RAN2并行化原理圖

當(dāng)需要產(chǎn)生random_num個(gè)隨機(jī)數(shù)并且可分配的線程數(shù)為thread_num時(shí)（thread_num＜N)，首先，通過2個(gè)初始種子計(jì)算出各段的初始值，然后各個(gè)線程將各段初始值作為種子開始依次計(jì)算random_num/thread_num個(gè)隨機(jī)數(shù)并輸出，直至產(chǎn)生夠random_num個(gè)隨機(jī)數(shù)為止.其并行化原理圖如圖1所示：

在各線程開始產(chǎn)生隨機(jī)數(shù)之前，需要根據(jù)2個(gè)初始種子先計(jì)算出每個(gè)段的初始值，即各線程的初始種子，由公式1.1可推導(dǎo)出:

如果將所有的線程從0號開始編程，并且0號線程的初始值即為初始種子seed0=（X0，Y0)，那么1號線程的初始值seed1=（d1×X0，d2×Y0)的計(jì)算方法如下所示：

其中L1，L2，…，L64為對應(yīng)于L的64位二進(jìn)制表示法中的每一位，取0或1.由公式1.3可推導(dǎo)得:

這樣，各個(gè)線程的初始值就都可以依次推出了，各個(gè)線程根據(jù)初始種子依次計(jì)算出自己所分配的隨機(jī)數(shù)序列，最終將各個(gè)線程產(chǎn)生的隨機(jī)數(shù)序列合并到一起就是所需產(chǎn)生的最終的隨機(jī)數(shù)序列.由于RAN2的周期約為2.3×1018，本文取L=246，即每個(gè)線程最多可以產(chǎn)生246個(gè)隨機(jī)數(shù)，并且該算法最多可以支持32684個(gè)線程.其并行算法的程序流程圖如圖2所示:

圖2 RAN2并行算法流程圖

其中，各個(gè)變量的含義為:

1)threads_num:線程數(shù);

2)random_num:總共要生成的隨機(jī)數(shù)的數(shù)量;

3)num:每個(gè)線程要生成的隨機(jī)數(shù)的數(shù)量;

4)id:線程的編號;

5)a[1]:Knuth38、RAN2、CombLec88的參數(shù)a1

6)a[2]:Knuth38、RAN2、CombLec88的參數(shù)a2

7)m[0]:Knuth38、RAN2、CombLec88的參數(shù)m1

8)m[1]:Knuth38、RAN2、CombLec88的參數(shù)m2

9)seed[i]i=0，1，2，3：初始種子，也是0號線程的初始種子

10)seed[id][0]和seed[id][1]:編號為id的線程的初始種子;

11)cout:移位計(jì)數(shù)器

12)seed[i][id] i=0，1:第id號線程的2個(gè)初始種子

13)result[random_num]:最終生成的隨機(jī)數(shù)序列.

由上圖可以得出，當(dāng)產(chǎn)生大量隨機(jī)數(shù)時(shí)，由于計(jì)算各個(gè)線程初始值和分配各線程任務(wù)時(shí)間很短故可忽略不計(jì)，假設(shè)各個(gè)隨機(jī)數(shù)發(fā)生器運(yùn)行一次的時(shí)間為單位1，則理論上當(dāng)采用thread_num個(gè)線程產(chǎn)生random_num個(gè)隨機(jī)數(shù)的時(shí)間約為O（random_num/thread_num），即相對于串行算法，其加速比應(yīng)為thread_num，但是由于線程的分配及共享變量的訪問等其他因素會導(dǎo)致加速比下降.

1.3基于MIC的并行化RAN2實(shí)現(xiàn)

Intel MIC（Many Integrated Core，集成眾核）架構(gòu)是英特爾公司專為高性能計(jì)算設(shè)計(jì)的、基于英特爾至強(qiáng)處理器和英特爾集眾核的下一代平臺[11].MIC在單一芯片上集成了約60個(gè)核，每個(gè)芯片計(jì)算峰值可達(dá)到每秒一萬億次以上的雙精度浮點(diǎn)運(yùn)算，處理復(fù)雜的并行應(yīng)用是MIC眾核架構(gòu)的優(yōu)勢.

MIC的應(yīng)用模式包括CPU原生模式、CPU為主MIC為輔模式、CPU與MIC對等模式、MIC為主CPU為輔模式和MIC原生模式五種.本文采用MIC原生模式進(jìn)行并行化，基于MIC的RAN2并行化實(shí)現(xiàn)分為三個(gè)步驟:首先，使用OpenMP[12，13]將串行程序并行化，通過TestU01測試保證并行化實(shí)現(xiàn)過程的正確性；其次，在CPU上進(jìn)行測試性能；最后，移植到MIC上做編譯選項(xiàng)優(yōu)化并進(jìn)行性能測試.

2　實(shí)驗(yàn)數(shù)據(jù)及性能分析

2.1 TestU01測試分析

目前常用的隨機(jī)數(shù)發(fā)生器測試庫有George Marsaglia 的Diehard Tests[14]，Donald Knuth的Classical Tests，美國國家科技標(biāo)準(zhǔn)局的NIST以及Pierre L`Ecuyer的TestU01[15，16].因TestU01包含共計(jì)216種455項(xiàng)測試，并且分為SmallCrush、Crush和BigCrush 3種測試方法，其測試范圍廣泛而且全面，因此本文將采用TestU01作為隨機(jī)數(shù)發(fā)生器測試庫.

本實(shí)驗(yàn)利用TestU01隨機(jī)數(shù)發(fā)生器測試庫分別對RAN2隨機(jī)數(shù)發(fā)生器的串行程序和4線程并行程序進(jìn)行了SmallCrush、Crush和BigCrush測試.圖3所示即為兩者的P-value統(tǒng)計(jì)分布情況，其中橫軸表示P-value的區(qū)間，縱軸表示測試結(jié)果處于不同區(qū)間的比例情況，.兩者均完全通過了TestU01中455項(xiàng)嚴(yán)格的測試，本文將串行和并行的P-value進(jìn)行了雙樣本KS統(tǒng)計(jì)測試，結(jié)果為0.4985，說明兩者服從同一分布.該結(jié)果在一定程度上反應(yīng)了RAN2隨機(jī)數(shù)發(fā)生器良好的隨機(jī)性并且證明了并行后隨機(jī)數(shù)發(fā)生器的正確性.

圖3串行與4線程并行程序TestU01測試P-value統(tǒng)計(jì)

2.2 CPU平臺下測試分析

本實(shí)驗(yàn)所選用的CPU平臺為兩個(gè)8核處理器Intel Xeon E5-2680 @ 2.7GHz，最多支持16個(gè)線程，在測試中分別對1、2、4、8和16個(gè)線程下產(chǎn)生1，000，000、10，000，000、100，000，000、1，000，000，000及10，000，000，000個(gè)隨機(jī)數(shù)的運(yùn)行時(shí)間進(jìn)行了測試，測試結(jié)果如表1所示:

表1基于CPU的時(shí)間測試

通過對表1時(shí)間測試的結(jié)果進(jìn)行分析，其加速比情況如圖4所示.

測試結(jié)果表明，當(dāng)所需產(chǎn)生的隨機(jī)數(shù)較少時(shí)，加速比并不明顯，這是因?yàn)楫?dāng)生成的隨機(jī)數(shù)較少時(shí)，總的運(yùn)行時(shí)間太短以致于初始化程序時(shí)間占較大比例，而各個(gè)線程并行產(chǎn)生隨機(jī)數(shù)的時(shí)間占用比例較小.因此隨著線程數(shù)的增加其加速比甚至有下降的情況；當(dāng)產(chǎn)生隨機(jī)數(shù)較多時(shí)，加速比并沒有隨線程數(shù)的增加而線性增長，這是因?yàn)榫€程數(shù)增加時(shí)，各線程和內(nèi)存之間交換的數(shù)據(jù)量增大，創(chuàng)建和銷毀線程將占用一部分系統(tǒng)資源及計(jì)算機(jī)系統(tǒng)對線程的調(diào)度管理都將降低并行效率.從圖中還可以看出，當(dāng)線程數(shù)相同時(shí)，隨著產(chǎn)生隨機(jī)數(shù)數(shù)量的增加，并行效率也隨之提高.例如，16個(gè)線程產(chǎn)生1，000，000個(gè)隨機(jī)數(shù)的加速比為1.547.產(chǎn)生10，000，000個(gè)隨機(jī)數(shù)的加速比為7.828，產(chǎn)生100，000，000個(gè)隨機(jī)數(shù)的加速比為10.837，產(chǎn)生1，000，000，000個(gè)隨機(jī)數(shù)的加速比為12.440，產(chǎn)生10，000，000，000個(gè)隨機(jī)數(shù)的加速比為14.273，因此該并行算法是可擴(kuò)展的.

圖4 CPU平臺下的加速比

2.3 MIC平臺下的測試及分析

本文所使用的MIC平臺為Intel Xeon Phi 3110P 57cores @ 1.10 GHz，首先在224個(gè)線程下產(chǎn)生1，000，000，000個(gè)隨機(jī)數(shù)，然后對各項(xiàng)編譯選項(xiàng)進(jìn)行時(shí)間測試，其中只有添加-O3（選擇最高的優(yōu)化級別“3”)編譯選項(xiàng)后優(yōu)化效果比較明顯，其相對-O2編譯選項(xiàng)的加速比為1.46，因此本文采用-O3編譯選項(xiàng).

本實(shí)驗(yàn)分別對1、2、4、8、16、32、56、112、168和224個(gè)線程下產(chǎn)生1，000，000、10，000，000、100，000，000、1，000，000，000 及10，000，000，000個(gè)隨機(jī)數(shù)的時(shí)間進(jìn)行了測試，測試結(jié)果如表2所示:

表2基于MIC的時(shí)間測試

通過對表2時(shí)間測試的結(jié)果進(jìn)行分析，可以獲得其加速比情況如圖5所示:

圖5 MIC平臺下的加速比

由圖6可知，基于MIC的RAN2并行化加速比效果較明顯，最優(yōu)加速比達(dá)117.07倍.當(dāng)產(chǎn)生隨機(jī)數(shù)較少時(shí)，各個(gè)線程并行產(chǎn)生隨機(jī)數(shù)的時(shí)間占總時(shí)間的比例較小，因此隨著線程數(shù)的增加并行效果并不明顯.當(dāng)產(chǎn)生隨機(jī)數(shù)越來越多時(shí)，性能提升效果也越來越明顯，但是加速比并沒有隨線程數(shù)的增加而線性增長，這是因?yàn)榫€程數(shù)增加時(shí)，各線程和內(nèi)存之間交換的數(shù)據(jù)量也隨之增大.此外，在MIC卡上設(shè)置的線程數(shù)并不是越多越好，線程數(shù)越多則將導(dǎo)致開銷越大.因此要根據(jù)所需產(chǎn)生隨機(jī)數(shù)個(gè)數(shù)設(shè)置合適的線程數(shù)以確保獲得最佳的加速比.例如圖6中基于MIC產(chǎn)生100，000，000個(gè)隨機(jī)數(shù)時(shí)，在56個(gè)線程左右的加速比最高，多于56線程時(shí)加速比出現(xiàn)下滑趨勢.此外，當(dāng)線程數(shù)相同時(shí)，隨著產(chǎn)生隨機(jī)數(shù)數(shù)量的增加，加速比逐漸呈增長狀態(tài)，因此該并行算法在MIC平臺上也是可擴(kuò)展的.

3　結(jié)束語

本文在分析了RAN2隨機(jī)數(shù)發(fā)生器串行算法的基礎(chǔ)上，設(shè)計(jì)并實(shí)現(xiàn)了在MIC平臺上的RAN2隨機(jī)數(shù)發(fā)生器的并行化方案，并且通過TestU01測試驗(yàn)證了其正確性，基于CPU平臺的性能測試和基于MIC平臺的性能測試結(jié)果表明該方法可以有效的減少RAN2的運(yùn)行時(shí)間，提高其產(chǎn)生隨機(jī)數(shù)的效率，使得RAN2隨機(jī)數(shù)發(fā)生器能夠更為廣泛的應(yīng)用于高性能領(lǐng)域的科學(xué)研究中.在下一步的工作中將進(jìn)行其它隨機(jī)數(shù)發(fā)生器的并行化研究.

參考文獻(xiàn)：

〔1〕Knuth D E. The Art of Computer Programming[M]. 2nded. Addis on-Wesley Publishing Company.2002：1-177.

〔2〕Knuth D E.計(jì)算機(jī)程序設(shè)計(jì)藝術(shù)(第2卷)：半數(shù)值算法[M].蘇運(yùn)霖，譯.3版.國防工業(yè)出版社，2002.8-163.

〔3〕D H LEHMER. Mathematical methods in large-scale

computing units [J]. Annals of the Computation Laboratory of Harvard University，1951：141-146.

〔4〕P.L'Ecuyer. Combined multiple recursive random number generators [J].Operations Research. Vol.44，No.5，Sep-Oct 1996：816-82.

〔5〕H.C.Tang. Combined random number generator via the generalized Chinese remainder theorem [J]，Journal of Computational and Applied Mathematics，Volume 142，Issue 2，15 May 2002：377-388.

〔6〕Press，W. H. and Teukolsky，S. A. 1992. Numerical Recipes in C.Cambridge University Press，Cambridge.

〔7〕楊尚琴.多層次并行算法與MPI－2新特性的研究及應(yīng)用[D].成都理工大學(xué)學(xué)位論文，2009.

〔8〕魏公毅，楊自強(qiáng).關(guān)于并行隨機(jī)數(shù)發(fā)生器的若干算法[J].數(shù)值計(jì)算與計(jì)算機(jī)應(yīng)用，2001(4)：311-320.

〔9〕MICHAEL MASCAGNI，ASHOK SRINIVASAN. Algorithm 806：SPRNG：a scalable library for pseudorandom number generation[J]. ACM Trans. Math. Softw，2000，26(3)：436-461.

〔10〕THOMAS BRADLEY，JACQUES TOIT，MIKE GILES，et al. Parallelisation techniques for random number ge nerators [J]. GPU Gems：Emerald Edition，2011：231-24..

〔11〕Wang Endong，Zhang Qing，Shen Bo，et al. High-Performance Computing on the Intel Xeon Phi-How to Fully Exploit MIC Architectures [M]. China WaterPower Press. 2012：79-81.

〔12〕Shameem Akhter，Jason Roberts著多核程序設(shè)計(jì)技術(shù)－通過軟件多線程提升性能[M].李寶峰，富弘毅，李韜譯，譯電子工業(yè)出版社，2007.

〔13〕周偉明.多核計(jì)算與程序設(shè)計(jì)[M].華中科技大學(xué)出版社，2009.

〔14〕George Marsaglia. The Marsaglia Random Number CDROM including the Diehard Battery of Tests of Randomness [M/CD] http：//www.stat.fsu.edu/pub/ diehard/.

〔15〕P L’ECUYER and R SIMARD，TestU01：A C Library for Empirical Testing of Random Number Generators [J]. ACM Transactions on Mathematical Software，2007，33(4)：Article 22.

〔16〕McCullough，B. D. A Review of TESTU01[Z]，Appl. Econom，to appera，2006.

中圖分類號：TP311.52

文獻(xiàn)標(biāo)識碼：A

文章編號：1673-260X（2015）09-0022-04