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

蔡曉龍（洛陽師范學院，河南　洛陽　471022）

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

蔡曉龍
（洛陽師范學院，河南洛陽471022）

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

關鍵詞：隨機數發生器；RAN2；并行化；MIC

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

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

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

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

1.1 RAN2串行算法

RAN2隨機數發生器是由2個相互獨立的LCG隨機數發生器組合而成的，通過將2個LCG發生器的結果進行合并就可以得到最終的RAN2隨機數.RAN2隨機數發生器的遞推公式如下:

參數為:

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

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

1.2 RAN2并行化實現

圖1 RAN2并行化原理圖

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

在各線程開始產生隨機數之前，需要根據2個初始種子先計算出每個段的初始值，即各線程的初始種子，由公式1.1可推導出:

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

其中L1，L2，…，L64為對應于L的64位二進制表示法中的每一位，取0或1.由公式1.3可推導得:

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

圖2 RAN2并行算法流程圖

其中，各個變量的含義為:

1)threads_num:線程數;

2)random_num:總共要生成的隨機數的數量;

3)num:每個線程要生成的隨機數的數量;

4)id:線程的編號;

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

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

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

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

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

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

11)cout:移位計數器

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

13)result[random_num]:最終生成的隨機數序列.

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

1.3基于MIC的并行化RAN2實現

Intel MIC（Many Integrated Core，集成眾核）架構是英特爾公司專為高性能計算設計的、基于英特爾至強處理器和英特爾集眾核的下一代平臺[11].MIC在單一芯片上集成了約60個核，每個芯片計算峰值可達到每秒一萬億次以上的雙精度浮點運算，處理復雜的并行應用是MIC眾核架構的優勢.

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

2　實驗數據及性能分析

2.1 TestU01測試分析

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

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

圖3串行與4線程并行程序TestU01測試P-value統計

2.2 CPU平臺下測試分析

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

表1基于CPU的時間測試

通過對表1時間測試的結果進行分析，其加速比情況如圖4所示.

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

圖4 CPU平臺下的加速比

2.3 MIC平臺下的測試及分析

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

本實驗分別對1、2、4、8、16、32、56、112、168和224個線程下產生1，000，000、10，000，000、100，000，000、1，000，000，000 及10，000，000，000個隨機數的時間進行了測試，測試結果如表2所示:

表2基于MIC的時間測試

通過對表2時間測試的結果進行分析，可以獲得其加速比情況如圖5所示:

圖5 MIC平臺下的加速比

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

3　結束語

本文在分析了RAN2隨機數發生器串行算法的基礎上，設計并實現了在MIC平臺上的RAN2隨機數發生器的并行化方案，并且通過TestU01測試驗證了其正確性，基于CPU平臺的性能測試和基于MIC平臺的性能測試結果表明該方法可以有效的減少RAN2的運行時間，提高其產生隨機數的效率，使得RAN2隨機數發生器能夠更為廣泛的應用于高性能領域的科學研究中.在下一步的工作中將進行其它隨機數發生器的并行化研究.

參考文獻：

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

〔2〕Knuth D E.計算機程序設計藝術(第2卷)：半數值算法[M].蘇運霖，譯.3版.國防工業出版社，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新特性的研究及應用[D].成都理工大學學位論文，2009.

〔8〕魏公毅，楊自強.關于并行隨機數發生器的若干算法[J].數值計算與計算機應用，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著多核程序設計技術－通過軟件多線程提升性能[M].李寶峰，富弘毅，李韜譯，譯電子工業出版社，2007.

〔13〕周偉明.多核計算與程序設計[M].華中科技大學出版社，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

文獻標識碼：A

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