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

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

Screen Shot 2013-11-14 at 2.29.39

 

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

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

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

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

“CUDA简单科学运算尝试:Make It Quick and Clean!”的一个回复

评论已关闭。