CUDA：蒙特卡洛以及随机数支持

由于随机数是整个计算物理的基石，实在太过于重要，干脆用一片文章来叙述CUDA有关于随机数产生器的种种。

目录

这一篇文章暂时更新一些随机数相关的东西，第一步的简单应用请参考我的另一篇文章，对于使用CUDA做蒙特卡洛积分的一些探讨

首先：学习资料

http://docs.nvidia.com/cuda/curand/introduction.html

文章中很多地方我就摘抄guide的原话了，避免因为翻译造成的理解误差。

作为一个严谨的编程者，很多事情不仅要考虑系统API的速度，也要对其严谨性审查，尤其蒙特卡洛这种靠随机数起家的玩意。好在CUDA给我们提供了一套非常靠谱的随机数产生器

The CURAND library provides facilities that focus on the simple and efficient generation of high-quality pseudorandom and quasirandom numbers.

–摘自introduction

这个quasirandom说的口气很大嘛！丁泽军老师的计算物理课也就敢说伪随机而已。

那我们来一探究竟吧！

首先，CUDA的随机数提供分为两种类型，第一是由GPU运算并且交给主机的随机，让我一下子就想到了计算物理课上的那个“随机数硬件产生器”，大抵就是如此了。

用过中国银行网银的孩子估计见过那种高大上的密码生成器，它采用每分钟给出一个6位数字的序列，并和网络服务器同步，有效期有几年，算下来这个给出随机数的小玩意要在其一生中把6位数的一半作为排序给出，并且要和网络服务器达到分钟的精度，那么必然是存储好的了。且不论其安全性（以后有机会可以研究下），但就这样一个超大规模又涉及金融安全的玩意，如果你拿16807算法来生成随机序列还真不怎么靠谱。

另一种就是在device上直接运行的随机数产生器，由于临近考试在敢作业，就先从device的随机数产生器讲起吧。

关于device api，NVIDIA这么叙述

To use the device API, include the file curand_kernel.h in files that define kernels that use CURAND device functions.The device API includes functions pseudorandom generation for and quasirandom generation .

CUDA采用几种随机数生成算法，包括：

 XORWOW， MRG32k3a， MTGP32

具体的算法Google吧。以后再分析。（谁刚刚还在说自己是严谨的编程者呢。。口亨。。。）

于是，随机数产生器一般要做且只需要两件事情，初始化与获得值，大约还记得Pascal的随机数产生器是这么写的

1. …
2. randomize();
3. a:=random();
4. …

非常省时省心，可对于CUDA我们就得一步步来咯。

初始化函数的定义如下：

1. __device__ void
2. curand_init (
3.     unsigned long long seed, unsigned long long sequence,
4.     unsigned long long offset, curandState_t *state)

原文解释如下

The curand_init() function sets up an initial state allocated by the caller using the given seed, sequence number, and offset within the sequence. Different seeds are guaranteed to produce different starting states and different sequences. The same seed always produces the same state and the same sequence. The state set up will be the state after

$2$

67 ⋅ sequence + offset calls to curand() from the seed state.

CUDA可以产生若干随机序列，每一个随机序列产生有两个属性：第一是产生其的种子，也就是seed，另一是随机数产生的偏移量，因为我们知道随机数产生都是一个迭代过程，CUDA提供了在初始化时候就迭代若干次的算法来解决种子过少的问题。

传统而常见的种子提供做法是使用当前系统时间（微秒数），来作为种子给出，对于串行程序这简单而有效，但是显然对于并行程序，这并非明智的做法，因为我们在主机启动GPU内核的时候时间是一定的且只能用这个值作为种子，即使进行一些处理也未必可以达到要求，于是乎，CUDA提供了这样一个解决方案，对于非严谨的使用者反而更加使用。

对于CUDA的每个随机序列，使用一个State给出标志，然后在调用某一序列的下一个值的时候使用

1. __device__ unsigned int
2. curand (curandState_t *state);

即可。给出例子，完整代码请参见我的GitHub

对于内核代码

1. __global__ void random_gpu(long* C,long* time,curandState*state)
2. {
3.     long i = threadIdx.x;
4.     long seed=(*time)*(i+1);//因为所有给定时间一定，所以我们只能通过对时间进行简单处理
5.     int offset=0;//完全独立的序列，所以offset全部为零来节约时间
6.     curand_init (seed,i,offset,>state[i]);//设置第i个随机序列
7.     C[i]=curand(>state[i]);//获得第i个随机序列的随机值
8. }

对于主函数

1. int main()
2. {
3.     size_t size = N * sizeof(float);
4.  
5.     long* C=new long[N];
6.     long st=getCurrentTime();
7.     curandState *state;
8.     long *d_C;
9.     cudaMalloc(>state,sizeof(curandState)*N);//设立随机状态列
10.     cudaMalloc(>d_C, size);
11.     random_gpu<<<1,N>>>(d_C,getCurrentTimeForDev(),state);
12.     cudaMemcpy(C, d_C, size, cudaMemcpyDeviceToHost);
13.     long ed=getCurrentTime();
14.     printf(“gpu running time:%ld\n”,ed–st);
15.     cudaFree(d_C);
16.     for(int i=0;i<10;i++)
17.     {
18.         printf(“%ld “,C[i]);
19.     }
20.     delete[] C;
21.     printf(“\n”);
22. }

于是就完成了，还有几个细节需要处理，因为CUDA略微繁琐的内存管理，所有输入值仅仅提供指针的形式给出，所以我们用getCurrentTimeForDev()来直接给出device指针的时间

1. long getCurrentTime()  
2. {  
3.     struct timeval tv;  
4.     gettimeofday(>tv,NULL);  
5.     return tv.tv_sec * 1000000 + tv.tv_usec;  
6. }
7. long*getCurrentTimeForDev()
8. {    long *time;
9.     cudaMalloc(>time,sizeof(long));
10.     long *timenow=new long;
11.     *timenow=getCurrentTime();
12.     cudaMemcpy(time,timenow,sizeof(long),cudaMemcpyHostToDevice);
13.     return time;
14. }

然后运行下：

1. CUDA_Learning$ ./rand
2. gpu running time:490969
3. 3500780692 144000514 797145604 4199381721 2915119786 2557417372 4151631408 1974695633 2578972200 2865224429 

随机数生成这事情确实慢的惊人。。。居然跑了半秒。

下面的样例直接切到用这个随机数序列来做蒙特卡洛积分了。请戳相应链接。