CUDA N-Body模拟

这学期的多核课简直爽,之前一直想玩的cuda也玩上了。真的是相见恨晚,我见识到什么叫吊。

CUDA:

CUDA就是Nvidia GPU自带的计算平台,它是为了克服GPGPU的弱点而诞生的。我不知道ATI有什么类似的,但CUDA肯定是Nvidia的杀手锏。因为他计算能力确实强,而且编程非常灵活简单。唯一的不足是需要对硬件有所了解。当然了解并行算法也是非常重要的。

简单说下CUDA的基本原理。CUDA是牺牲单个核的主频,以核的数量来取胜。比如我笔记本的GT540M就有96个CUDA core。也就是说最多96个线程同时运行。基于这点CUDA的自然运用就是大量线程,而在CUDA里,线程的切换基本不消耗资源。问题在于如何组织这些线程。可以看basics of cuda。CUDA每个任务称为一个kernel,一个kernel包含一个Grid,一个Grid可以有很多Block,Block中可以有很多线程。而且Grid和Block可以是3维的。考虑一维情况,一个kernel的总线程是GridDim*BlockDim。而CUDA的调度单元是一个warp,一个warp有32个线程。这的意思并不是32个线程同时执行一条指令。比如我有16个CUDA core,我是先执行前16个线程,然后后16个线程。而且最近CUDA还支持指令并行,比较复杂,详见simt-architecturestack overflow.

N-Body:

N-body说白了,就是万有引力。所以计算复杂度是O(N^2),(虽然有O(Nlog(N))的算法),1000个粒子,对CPU来说已经负担很大了。所以要使用GPU。所以算法很简单,关键是怎么优化。如果是最原始的方法,每一个线程中对都计算N个粒子对其作用力,那要读取N*N次内存。所以要采用分块的方式计算,P个线程协作,每次读取P个到共享内存,然后计算P个之间的作用力,那么只需N*N/P次内存读取(共享内存很快)。然后为了计算简单,这里没有碰撞,假定粒子能穿过粒子,其作用力为 sum(ri/(ri^2+eps^2)^(1.5)) ,eps为一个小量,防止作用力无穷大。

实现:

为了简单,画图用python的matplotlib,计算用CUDA。这里涉及接口问题,可以自己了解一下。最后的效果是模拟二维的。如下图的双星也是很漂亮(初始到一段时间后)。一共4096个粒子,30帧每秒,相当于(30*4096*4096*20 / 1024 / 1024 / 1024)  = 9GFLOPS,CPU顶不住。

1.two_cluster22.two_cluster3.two_cluster1