基于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)。宋博文，硕士生，主研领域：高性能计算。周津羽，研究员。华诚，工程师。刘逍，硕士生。周晓辉，教授。