基于MIC的MRG32k3a并行化設計與實現

宋博文　周津羽　華　誠　劉　逍,　周曉輝,

1(西安郵電大學計算機學院　陜西 西安 710121)

2(陜西省高性能計算研究中心　陜西 西安 710121)

?

基于MIC的MRG32k3a并行化設計與實現

宋博文1周津羽2華誠2劉逍1,2周曉輝1,2

1(西安郵電大學計算機學院陜西 西安 710121)

2(陜西省高性能計算研究中心陜西 西安 710121)

摘要隨機數產生器在工程模擬等領域獲得廣泛應用，MRG32k3a是一種性能優異的隨機數產生器，但產生速率較慢。針對這種情況，在研究MRG32k3a串行算法的基礎上，利用算法并行化理論，提出一種基于MIC(Many Integrated Core)平臺的MRG32k3a并行化方法。實驗結果表明，該方法能通過TestU01的全部測試，移植到MIC平臺后加速比與線程數呈線性增長關系，相對CPU單線程的最佳加速比為17.73。

關鍵詞隨機數產生器MIC并行化MRG32k3aTestU01

DESIGN AND IMPLEMENTATION OF MRG32K3A PARALLELISATION BASED ON MIC

Song Bowen1Zhou Jinyu2Hua Cheng2Liu Xiao1,2Zhou Xiaohui1,2

1(School of Computer Science and Technology,Xi’an University of Posts and Telecommunications,Xi’an 710121,Shaanxi,China)2(Shaanxi Research Center for High Performance Computing,Xi’an 710121,Shaanxi,China)

AbstractRandom number generator is widely used in engineering simulation, MRG32k3a is a random number generator with excellent performance, but its generation rate is slow. In view of this, this article presents an MIC platform-based MRG32k3a parallelisation approach using the theory of algorithm parallelisation according to the study on MRG32k3a serial algorithm. Experimental results show that the approach can pass all the tests of uniform random number test-library TestU01. After transplanting to MIC platform, the relationship of the speedup ratio and the number of threads increases linearly, and the best speedup ratio relative to CPU single-thread reaches 17.73.

KeywordsRandom number generatorMany integrated core (MIC)ParallelisationMRG32k3aTestU01

0引言

隨機數產生器是用來產生隨機數的裝置[1]，一般分為真隨機數產生器和偽隨機數產生器。真隨機數通過物理實驗或自然噪音等方式產生，成本高；偽隨機數通過若干種子利用數學算法遞推產生周期性的隨機數序列，不需要外部物理硬件的支持，并且具有類似于真隨機數的統計特征；因此，在科學研究和工程模擬領域主要應用偽隨機數產生器。相對于現在普遍應用的偽隨機數產生器LCG(Linear Congruential Generator)[2]，MRG32k3a[3-5]具有長周期、隨機數序列質量優異的特點，但較慢地產生速率制約了它的應用。

隨著計算機系統結構和計算機網絡技術的快速發展，科學研究和工程模擬對隨機數產生器的性能及速率要求也日益增長，研究隨機數產生器并行化是迫切的、必然的。并行計算技術迅速發展[6]，將其應用于隨機數產生器，成為提高隨機數產生器產生效率的方法之一。目前國內外已經做了一些隨機數產生器的并行技術的研究，并行方法總體分為兩種，第一種是每個線程獨自應用不同類別的隨機數產生器，第二種是將一個隨機數產生器并行化到每一個線程中。第一種方法易于并行化，2000年Michael Mascagni等研發出了基于第一種方法的并行化隨機數產生器庫SPRNG(Scalable Parallel Random Number Generators Library)，每個線程所使用的產生器都不相同[7]；后來Shuang Gao等在此基礎上研發出基于統一計算設備架構CUDA(Compute Unified Device Architecture)平臺的并行化隨機數產生器庫GASPRNG(GPU accelerated scalable parallel random number generator library)[8]；Kinga Marton等于2011年提出的一種并行化隨機數產生器，各線程混合使用了不同類型的隨機數產生器[9]。第二種并行化方法則較難實現，優點是能保持隨機數產生器的原始序列。 2011年，Thomas Bradley等將第二種并行化方法總結為Simple skip ahead、Strided skip ahead和Hybrid三類，并且描述了基于GPU的MRG32k3a、Sobol和MT19937三種隨機數產生器的并行化算法及性能分析[10]。L’Ecuyer等人2001年也早以概述了MRG類隨機數產生器的并行化的理論過程，但沒有驗證并行化后的結果正確性和運行效果[11]。

目前隨機數產生器的并行化研究工作主要集中于多核中央處理器CPU平臺，缺少基于Intel最新產品MIC平臺[12]并行化的相關理論依據和性能分析。另外，以上文獻均沒有詳細闡述MRG32k3a并行化的具體使用方法、并行化實施方案及并行化結果分析。本文利用“分而治之”思想實現MRG32k3a并行化，設計了一種簡單的并行化規則，并確定了相應的參數，實現了按周期“分而治之”并行產生隨機數的算法，更主要的工作是基于MIC進行了MRG32k3a并行化后的性能測試分析，檢驗其可行性和優越性。文獻[10]中Thomas Bradley總結的Simple skip ahead方法實質上就是 “分而治之”思想在隨機數產生器應用上的體現，在一個周期內劃分線程和分配任務，每個線程獨立產生一段周期內的連續的隨機數子序列。

本文主要研究工作分為以下三個階段：1) 理論研究MRG32k3a的串行與并行算法，確定最有效的并行化方案，進行理論加速比和可擴展性分析；2) 在Intel Xeon E5-2680上進行了MRG32k3a串行和并行的大量開發與測試，其中選用最為完備的均勻隨機數產生器測試庫TestU01進行隨機數序列性能測試[13,14]，借助TestU01檢測保證了并行化過程的正確性，最終MRG32k3a串行及并行算法均能完全通過TestU01 測試，同時完成相應時間測試和加速比分析；3) 移植到MIC平臺，通過測試挑選出最優的Intel編譯選項，然后完成時間測試和性能分析。實驗數據顯示MRG32k3a并行化后的加速比和并行效率理想，MIC相對CPU單線程的最佳加速比為17.73。

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

1.1MRG32k3a串行算法

組合多重遞歸隨機數產生器[3,4]CMRG( Combined Multiple Recursive Generator)遞推公式為：

xj,n=(aj,1xj,n-1+ … + aj,kxj,n - k)modmj

(1)

為兼顧速度和長周期性,j= 2、k= 3時比較理想[5]。MRG32k3a屬于此種類型的CMRG，遞推公式為：

(2)

其中n≥3，a1,2=1 403 580,a1,3=-810 728,a2,1=527 612,a2,3=-1 370 589,m1=232-209,m2=232-22 853。

產生[0,1)之間的均勻隨機數un。

zn= (x1,n+ x2,n)modm1

(3)

1.2MRG32k3a并行化算法設計

MRG32k3a產生的隨機數序列之間存在相互依賴關系，最大維度只有45[11]，采用分治法可以避免維度約束。每個線程跳到隨機數序列特定位置后連續產生隨機數，意思是說把原始隨機數序列分割為V段，每個子序列都是原始序列的連續子序列，這些特定位置是每個線程初始值。本文將闡述MRG32k3a利用分治法并行化的具體實現過程。

設隨機數序列每個狀態為向量對形式：

Xi,n=(xi,nxi,n +1xi,n + 2)T

MRG32k3a遞推公式則轉化為：

Xi,n +1=AiXi,nmodmii=1,2

(4)

狀態向量組Xi,n+w可以直接由Xi,n計算：

(5)

令p是MRG32k3a的周期，T為轉換函數，則T(Xn)=Xn+1，T(Xp)=X。整個周期的隨機數序列劃分成V=2v個長度為W=2w的隨機數子序列，其中v+w=191。X0是隨機數序列的初始狀態，每個隨機數子序列的初始狀態定義為：

C1=X0,C2=Tw(X0),…,Cg=Tw(Xg-1)=T(g-1)w(X0)

本文選取w=127,v=64，即最多運行V=264個線程，各線程最多產生W=2127個隨機數序列，圖1為相應的并行化原理圖。

圖1　MRG32k3a并行化原理圖

如圖1所示，MRG32k3a并行化的關鍵點是計算每個線程的初始狀態Cg。其中，三個參量乘積的模運算[1]定義為：

X·Y·Zmodm=((X·Ymodm)·Z)modm

(6)

MRG32k3a并行化的具體算法描述如算法1：

算法1MRG32k3a并行化算法

輸入： 種子seed[6]，線程數thread_num，任務量random_num

輸出： random_num個隨機數u

Begin

(1) for i = 0 to thread_num do

(1.1) for j = 0 to 6 do

C[i][j] = seed[j]

end for

end for

(2) for all id where 0≤id≤thread_num do

(2.1) for j = 0 to random_num/thread_num do

p1 = a12 * C[id][1] - a13 * C[id][0]

p2 = a21 * C[id][5] - a23 * C[id][3]

u = ((p1 > p2) ?

(p1 -p2)/m1 : (p1 -p2 + m1)/m1

end for

end for

End

MRG32k3a串行算法產生n個隨機數序列，所需時間與n大小成正比，因此時間復雜度為O(n)或線性時間，它的并行算法在g個線程情況下時間復雜度為O(n/g)，所以其絕對加速比為O(g)，有效利用了線程數量，加速比隨著線程數增加而線性增長。MRG32k3a串行算法的任務量與產生的隨機數總量random_num成正比，同時MRG32k3a的并行算法是將總任務量平均分配給g個線程，理想情況下，并行算法的加速比接近于線程數。而當任務量較少時，一定的通信開銷會額外消耗時間，比如分配線程、訪問共享變量、內存間轉換等操作，相應的加速比會下降，而任務量逐漸增加后，通信開銷所占比例越來越小，逐漸呈現出線性加速。

1.3基于MIC并行算法實現

Intel推出的基于集成眾核MIC架構的至強融核系列產品是為了解決高度并行計算問題。它基于x86架構，支持OpenMP、pThread等并行編程模型。英特爾MIC架構具有更小的內核和更多的硬件線程，以及更寬的矢量單元，是提高整體性能、滿足高度并行化應用需求的理想之選[12]。

并行程序執行時，多線程之間相互通信會使程序運行變得復雜，選擇合適的并行粒度才能最大限度地提高并行的效率。MRG32k3a并行化屬于細粒度并行，每線程單獨產生一段連續隨機數子序列，并發度高、占用資源較少、代碼較短、各線程間數據交換少，適合移植到MIC上運行。

本文基于MIC的MRG3ak3a并行化實現分為三個步驟：首先，將串行程序使用OpenMP并行化，通過TestU01測試其隨機數序列的準確性；再次，測試CPU上并行版本的性能，如果加速比不能隨著核數增加而線性增長，則需要對程序進行優化，使之發揮出多核的威力；最后，移植到MIC上進行編譯選項篩選，Intel的各種優化編譯選項在不同程度上會提升運行效率，為最大發揮MIC性能，要選擇出最適合MRG32k3a并行程序運行的編譯選項，接著完成MIC上的性能測試。

2實驗數據及性能分析

2.1TestU01測試分析

為保證并行化結果的正確性，本文選用均勻隨機數統計檢測庫TestU01進行隨機數性能測試。TestU01相比Diehard，NIST等統計檢測庫，其檢測最為嚴格，涉及216種、455個統計檢測。綜合SmallCrush、Crush、BigCrush的三大類測試結果，說明了MRG32k3a串行與并行程序能夠全部通過TestU01測試。串行版本與4線程并行版本的TestU01測試的P-value統計分布情況，如圖2所示。圖2中串行版本P-value落在區間(0.08，0.92)之間比例為88.56%，4線程并行版本P-value落在該區間的比例為87.24%，通過的檢測比例基本相同。對串行和并行程序的測試數據進行了雙樣本柯爾莫諾夫-斯米爾諾夫KS(Kolmogorov-Smirnov)統計測試，P-value為0.6517，進一步說明它們測試數據近似服從相同分布。

圖2　MRG32k3a串行與4線程并行程序TestU01測試的P-value統計分布

2.2基于CPU的并行化性能測試

實驗環境：搭建的CPU平臺為兩個8核處理器Intel Xeon E5-2680 @ 2.7 GHz。

實驗數據及分析：為更好地評估MRG32k3a程序的并行化成果，首先基于CPU進行了基準測試，在1、2、4、8、16線程情況下產生不同數量隨機數的時間測試結果如表1所示。其中產生同等數量隨機數的條件下，若線程數成倍增加，測試時間則成倍縮短；同等線程數條件下，若產生隨機數的數量成倍增加，測試時間也成倍增長；而且一定情況下，產生隨機數量越多其可擴展性越好。圖3表示相應線程下的加速比，橫軸為線程數，縱軸為加速比，基于CPU的MRG32k3a并行化加速比與線程數呈線性增長，此時具有一定擴展性，16線程時達到12.1609倍的加速比，而且在2、4、8線程時并行化效率接近100%，因此MRG32k3a并行程序可以移植到MIC。

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

續表1

圖3　基于CPU的加速比趨勢圖

2.3基于MIC的并行化性能測試

實驗環境：本文MIC平臺使用的是單個MIC卡Intel Xeon Phi 3110P 57cores @ 1.10 GHz，它含有57個 x86 架構指令集的核心，每個物理核心帶有4線程。

篩選優化選項：首先在224線程產生1010個隨機數情況下，依次添加不同編譯選項進行時間測試，測試結果如表2所示，圖4為相應的加速比。

表2　不同Intel編譯選項的時間測試(s)

圖4　不同編譯選項的時間測試加速比

如圖4所示，基于MIC并行化添加編譯選項-O3、-no-prec-div 、-fimf-domain-exclusion時加速比有所上升，優化效果比較明顯。 其中-O3表明選用最高的優化級別“3”，-no-prec-div指使用乘倒數替代除法來簡化運算，-fimf-domain-exclusion指排除與算法無關的浮點運算，減少非正常浮點操作對性能的影響。

實驗數據及分析：篩選出合適的編譯選項后，基于MIC平臺實現1、2、4、8、16、32、56、112、168、224線程產生不同數量隨機數序列的時間測試，測試結果如表3所示，圖5為相應線程下的加速比。

表3　基于MIC的時間測試(s)

圖5　基于MIC的加速比趨勢圖

如圖5所示，基于MIC的MRG32k3a并行化加速比效果也非常明顯，最大達到216.2105，部分線程的并行效率超過了100%。當產生隨機數越多時，提升效率越明顯。另外，在MIC卡上設置的線程數并不是越多越好，線程數太多，開銷會比較大。因此要設置合適的線程數以確保MIC核的高利用率，比如圖5中基于MIC產生107個隨機數時在56線程左右的效率最好，多于56線程時加速比出現下滑趨勢。當產生隨機數數量不斷增加，如產生1010數量的隨機數，加速比逐漸呈線性狀態，此時基于MIC也具有一定可擴展性。

由于C++語言標準庫Boost、隨機數產生器庫SPRNG中沒有隨機數產生器MRG32k3a，本文選取了英特爾數學核心函數庫Intel MKL(Intel Math Kernel Library)中MRG32k3a的并行化結果進行性能對比。如表4所示，測試數據來自文獻[10]，平臺為Intel Xeon E5410，pts/s指每秒可產生的隨機數的數量，本文中基于MIC實現MRG32k3a并行化的最優情況下每秒大約產生8.25E+08個隨機數，而Intel MKL 基于Intel Xeon E5410平臺4線程情況下每秒大約也只能產生1.04E+08個隨機數。另外，對比表1中基于兩個8核CPU產生隨機數的時間，單個MIC卡相對CPU單線程的最佳加速比為17.73，說明基于MIC實現MRG32k3a并行化具有一定優越性，雖然MIC單核的計算

能力弱于同期的CPU，但利用眾核并行的優勢，加速比依然理想。

表4　與Intel MKL中MRG32k3a的性能對比

3結語

本文主要研究了隨機數產生器MRG32k3a利用“分而治之”思想的并行化實現，將周期內隨機數序列分段產生，并行程序全部通過了TestU01的455個測試，并基于MIC進行時間測試和性能分析，最終加速比與線程數接近正比關系，部分線程并行效率超過100%，并具有一定擴展性。與串行算法相比，基于MIC的MRG3ak3a并行化后，運行效率得到極大改善，可廣泛應用于高性能領域的科學研究。我們下一步工作將集中于使MRG32k3a并行算法更加完善，嘗試CPU+MIC協同并行化，也會繼續研究討論其它隨機數產生器并行化實現。

參考文獻

[1] Knuth D E.The Art of Computer Programming[M].Volume 2:Seminumerical Algorithms,1981.

[2] L’Ecuyer P.Random numbers for simulation[J].ACM Transactions on Modeling and Computer Simulation,1990,33(10):85-97.

[3] L’Ecuyer P,Blouin F,Couture R.A search for good multiple recursive random number generators[J].ACM Transactions on Modeling and Computer Simulation,1993,3(2):87-98.

[4] L’Ecuyer P.Combined multiple recursive random number generators[J].Operations Research,1996a,44(5):816-822.

[5] L’Ecuyer P.Good parameters and implementations for combined multiple recursive random number generators[J].Operations Research,1999,47(1):159-164.

[6] 陳國良.并行算法的設計與分析[M].北京:高等教育出版社,2002.

[7] Mascagni M,Srinivasan A.Algorithm 806:SPRNG:A scalable library for pseudorandom number generation[J].ACM Transactions on Mathematical Software,2000.

[8] Gao S,Peterson G D.GASPRNG:GPU accelerated scalable parallel random number generator library[J].Computer Physics Communications,2013,184(4):1241-1249.

[9] Marton K,Suciu A,Petricean D.A parallel unpredictable random number generator[C]//Roedunet International Conference (RoEduNet),2011 10th.Piscataway:IEEE,2011:1-5.

[10] Bradley T,du Toit J,Giles M,et al.Parallelization techniques for random number generators[J].GPU Computing Gems,2010:231-246.

[11] L’Ecuyer P,Richard Simard E,Jack Chen,et al.An object-oriented random number package with many long streams and substreams[J].Operations Research,2001,50(6):1073-1075.

[12] 王恩東,張清,沈鉑,等.MIC高性能計算編程指南[M].北京:中國水利水電出版社,2012.

[13] L’Ecuyer P,Simard R.TestU01:A C Library for Empirical Testing of Random Number Generators[J].ACM Transactions on Mathematical Software,2007,33(22):1-40.

[14] L’Ecuyer P,Simard R.On the performance of birthday spacings tests for certain families of random number generators[J].Mathematics and Computers in Simulation,2001,55(1-3):131-137.

中圖分類號TP311.52

文獻標識碼A

DOI:10.3969/j.issn.1000-386x.2016.02.058

收稿日期：2014-06-02。陜西省自然科學基礎研究計劃項目(20 13JM8028)。宋博文，碩士生，主研領域：高性能計算。周津羽，研究員。華誠，工程師。劉逍，碩士生。周曉輝，教授。