CUDA简单科学运算尝试：Make It Quick and Clean!

使命召唤5中，当主角杀入丛林，寂静的只有月光，队长闻到了一丝不对劲的气息，于是告诉大兵们，

make it quick and clean!

这也差不多是并行运算编写的根本所在，快速，干净，同时一定不留1k脏内存。而且在多次的惨痛教训中浩浩发现，如果不把cuda的device程序写干净的话，你的电脑当然不会死机，cpu照常运转，只不过，屏幕别指望他再动了。

目录

计算物理作业拖了好久也改写了。。。还不知道怎么给丁老师解释呢，干脆做的漂亮点，用用cuda来算点小东西，废话不多说，我们就用cuda来算积分吧

首先基于Monte Carlo的积分运算完全是靠数量拼精度，虽说漂亮的抽样方法可以提高计算的精度，但是对于一个3维积分，如果你一维都踩不到1024个点你好意思谈精度？

不过可以算一下一维1024个点是什么概念 ，也就是$$2^{30}=1073741824$$个点。。然后再奢求下高精度。用CPU算，甭管你是mpi还是什么。。想要秒解的老OIer可以洗洗睡了。

于是，人生苦短，python太慢，我们用cuda(nvcc)。

对此，我归结出来一句话:

离效率越近，离灵魂越远。

但是为了计算，也只好把灵魂卖给指针他老人家了。

其实CUDA在很多方面和Monte Carlo简直是一对，漂亮的高度并行化，完美的模型，以及，高度的随机数支持。

Monte Carlo一维积分的算法就太简单了

随机取点取点取点。。。然后算算算。。。然后加起来平均平均平均。。。

那么，废话少说，CUDA走起，好，为了极大可能的优化效率我们先考虑这样一件事情:

假设我们才样N个点，期望的精度是k位，那么我们每个点的精度应该是多少。

根据平均性能考虑，在取样点$$k-Log_{10}N$$后面的位数实际上都是毫无意义的，也就是说，我们样点精度要求在$$k-Log_{10}N$$后面的都会为涛涛数据所淹没掉，对于double而言，精度在15-16位左右，对于float，精度在6-7为左右，而窝们的cuda代码打算用$$2^{30}$$来算积分（杀鸡焉用牛刀？不过我就用了），那么采样时使用float是科学。不过为了省事情。。。就干脆全部用double了

这是我们要首先处理的函数，$$\sqrt{x+\sqrt{x}}$$，写在c里面很简单：//不过后面会有些问题

1. double func0(double x)
2. {
3.     return sqrt(x+sqrt(x));
4. }

时间有限，废话少说，CUDA传入函数指针会很轻易造成系统崩溃（CPU还在运转），具体原因估计是入口只有一个之类的，鉴于我对函数指针了解并不丰厚，不在这里妄谈了。留个坑以后补上。

1. __global__ void inte_cell(double* l,double* r,double *res,long *time,curandState *state)
2. {
3.     //不要传函数指针
4.     int i = blockIdx.x*BlockN+threadIdx.x;
5.     long seed=(*time)+(i);//因为所有给定时间一定，所以我们只能通过对时间进行简单处理
6.     int offset=0;//完全独立的序列，所以offset全部为零来节约时间
7.     curand_init (seed,i,offset,>state[i]);//设置第i个随机序列
8.     double x=1,sum=0;
9.  double k=1;
10.     for(int j=0;j<k*Dev_Loop;j++)
11.     {
12.         x=(r[i]-l[i])*curand_uniform_double(>state[i]);
13.         sum+=sqrt(x+sqrt(x));//func(x);
14.     }
15.     res[i]=sum/Dev_Loop/k;
16.     __syncthreads();
17.  
18. }

恩，很简单嘛！就是一个函数，然后crazy采样加和即可。

关于随机数部分，请参见蒙特卡洛。

有几个小问题，基于安全性考虑，我们的l,r并没有公用同一空间。同时也是为下一种复杂采样的蒙特卡洛做准备。每个区间取值若干次求和平均。

在初步的测试中我才用了256*1024个计算细胞，这并不算是一个大的数据量，不过基于精细的为以后发展的考量，研究一种可用于G级别数据的加速加法还是有用的，于是我定义了这样一个玩意

1. __global__ void big_plus(double*a,double *res,int *threadNum)
2. {
3.     //为了尽可能的利用并行效率，加法采用两次树形相加的形式，每次加addNum个
4.     //如此可以对付2^30次的快速相加
5.     //虽然嘛。。。。这是毫无意义的啦！因为本程序只有2^20次的相加
6.     //不过留个接口以后用总是好事情
7.  
8.     int i=blockIdx.x*(*threadNum)+threadIdx.x;
9.     double sum=0;
10.     int k=i*addNum;
11.     for(int j=0;j<addNum;j++)
12.     {
13.         sum+=a[k];
14.         k++;
15.     }
16.     res[i]=sum/addNum;
17.     __syncthreads();
18.     
19. }

简单的说就是把大数据的相加也转化为一个高度并行的相加-合并-相加过程。

学过动归的孩子应该做过一道题：合并傻子，其实说来我当时做这个小函数的时候第一个想到的是这个….不过我们暂时先把这个放一边，因为这个函数的使用需要应用到CUDA中很重要的一个概念：流管理。我们现在先暂时不涉及这一话题。

更进一步的，我们需要几个小函数来处理一些细节问题，首先是关于初值给定的问题，这其实是满繁琐的事情，无非是给显存内的一个数组赋值为某一值…..就干脆写俩小函数解决。。。为了方便记忆就不重载了

1. double *DevValueD(double v,int len)//把host值转化为dev指针值
2. {
3.     double*res;
4.     cudaMalloc(>res,sizeof(double)*len);
5.     double *val=new double[len];
6.     for(int i=0;i<len;i++)
7.         val[i]=v;
8.     cudaMemcpy(res,val,sizeof(double)*len,cudaMemcpyHostToDevice);
9.     return res;
10. }
11. int *DevValueI(int v,int len)
12. {
13.     int*res;
14.     cudaMalloc(>res,sizeof(int)*len);
15.     int *val=new int[len];
16.     for(int i=0;i<len;i++)
17.         val[i]=v;
18.     cudaMemcpy(res,val,sizeof(int)*len,cudaMemcpyHostToDevice);
19.     return res;
20. }

于是一个基于蒙特卡洛的积分程序就很快搭建起来了（虽然足足折腾了两天，因为各种并发性的了解不足，典型代表是因为big_plus函数造成大量零的结果，或者是函数指针使得程序直接崩溃。）

1. int work()
2. {
3.     int threadPerBlock=BlockN;
4.     int numBlocks= 256;
5.     size_t size = BlockN *numBlocks*sizeof(double);
6.  
7.     long st=getCurrentTime();
8.  
9.     curandState *state;
10.     cudaMalloc(>state,sizeof(curandState)*1024*1024);//设立随机状态列
11.  
12.     double* d_A,*add_tem0,*add_tem1,*res;
13.     cudaMalloc(>d_A, size);
14.     cudaMalloc(>add_tem0, size/16);
15.     cudaMalloc(>add_tem1, size/256);
16.     cudaMalloc(>res,sizeof(double));
17.  
18.     inte_cell<<<numBlocks,threadPerBlock>>>(DevValueD(0.0,numBlocks*threadPerBlock),DevValueD(1.0,numBlocks*threadPerBlock),d_A,getCurrentTimeForDev(),state);
19.  
20.     double *result=new double[numBlocks*threadPerBlock];
21.     
22.     double fin_res=0;
23.     cudaMemcpy(result,d_A, size, cudaMemcpyDeviceToHost);
24.     for(int i=0;i<numBlocks*threadPerBlock;i++)
25.     {
26.         fin_res+=result[i];
27.     }
28.     fin_res/=(numBlocks*threadPerBlock);
29.     long ed=getCurrentTime();
30.  
31.     printf(“GPU running Time:%ld\n”,ed–st);
32.     printf(“final:%16.14f\n”,fin_res);
33.  
34.     cudaFree(d_A);
35.     cudaFree(d_A);
36.     cudaFree(add_tem0);
37.     cudaFree(add_tem1);
38.     cudaFree(state);
39.     cudaFree(res);
40. }
41.  

 

1. 更加详细的代码参见我的GitHub,里面包含对于CPU测试效率的代码。

还有就是关于CUDA的一些学习心得.

首先，任何没有__device__,__global__的函数都是在设备代码中无法调用的，而且，更进一步的是任何__global__的函数必须有3.5的计算等级才可以在设备代码中调用。

其次，函数指针，慎用。关于CUDA初步的体会是一切要非常的干净。

学过Haskell的可能可以体会的到那种严格到令人窒息的干净的编程（所以很多人都非常讨厌Haskell的io，因为她太违和了），在CUDA这种高度并行化的编程中，我们不得不把每个程序当做void格式的无副作用的纯函数，很奇怪吧。但是我们不得不把每个计算细胞的输入局限于输入的数组的一行，同样输出的一行。而且避免胡乱更改输入数据，在一大堆繁琐的接近硬件的代码中，获得一个漂亮的编程模型并不是一件容易的事情，所以有时候强迫症也是很重要的。

CUDA更像是为物理学家和工程师设计的。对于数学家以及合格的码农而言，使用递归思考就像呼吸一样自然，但CUDA中我并没有找到一套完美的，自动分配并行子程序的对于递归的解决方案，这或许会是我对于CUDA探索和自己创造的一个重点：

写一个基于CUDA的，完美的自并行递归模型。

当然CUDA的强大让我又回想起了另外一套我构想已久的玩意，

一种基于纯粹方程，自由度约束和强大通用概率求解内核的编程语言。

一套更像是物理而非数学的编程语言，靠蒙特卡洛，演化计算等框架为内核。当然得等我补完作业追到女神以后再茶余饭后写咯。

QUICK

首先必须明了的是CUDA无论模型多么抽象，GPU的核心数量总是有限的。启动更是需要时间的。所以我们需要一些测试已经经验来进行我们的编程工作。

上面写好的这个程序我们运行下，

1. 5–7$ ./integration 
2. GPU running Time:2282524
3. final:66.86583087598933
4. cpu:time:6758037,res:       0.942827

恩，速度还不错，彪CPU三倍。

这样，为了测试256*1024个线程的启动缓慢到底对程序影响有多大，我们对GPU核心的主循环乘一个系数，从0.5-16(2倍一个阶段)

1. 0.5 1726350
2. 1 2085732
3. 2 2656851
4. 4 3834121
5. 8 6236095

绘图，非常干净的线性向我们展示了出来

 

恩，前面的截距就是启动时间了（此处测试我们忽略了后面加和的时间）可见对于2*1024次加分系统启动等工作是非常耗费时间的，至于6-8则可以忽略了。

另外是当倍率比较高的时候，程序会不由自主的陷入崩溃（在Macbook上这是一件很恐怖的事情，我不得不反复强制关机），所以慎用。 

前面是简单的基于蒙特卡洛方法的直接积分，后面还会有些基于变换的。等我考完五个小时以后的理论力学吧。。。。。。

所有代码已更新到GitHub,请自行参阅。