目录
Step1：申请CPU和GPU内存空间并对数据进行初始化和拷贝操作。
N-body问题
原理
N-body问题（或者说N体问题），是一个常见的物理模拟问题。在N-body系统中，每个粒子体都会与剩下的其他粒子产生交互作用（交互作用因具体问题而异），从而产生相应的物理现象。
串行代码
首先，对每个点的位置和速度进行随机初始化。
|
1 |
void randomizeBodies(float *data, int n) { //初始化函数 for (int i = 0; i < n; i++) { data[i] = 2.0f * (rand() / (float)RAND_MAX) - 1.0f; } } |
然后，根据位置进行计算相互的作用力，以此为依据更新点的速度，该代码被整合到一个bodyForce函数中。
|
1 |
void bodyForce(Body *p, float dt, int n) { //计算相互作用力并记录速度 for (int i = 0; i < n; ++i) { float Fx = 0.0f; float Fy = 0.0f; float Fz = 0.0f; for (int j = 0; j < n; j++) { float dx = p[j].x - p[i].x; float dy = p[j].y - p[i].y; float dz = p[j].z - p[i].z; float distSqr = dx*dx + dy*dy + dz*dz + SOFTENING; float invDist = rsqrtf(distSqr); float invDist3 = invDist * invDist * invDist; Fx += dx * invDist3; Fy += dy * invDist3; Fz += dz * invDist3; } p[i].vx += dt*Fx; p[i].vy += dt*Fy; p[i].vz += dt*Fz; } } |
根据新的速度，更新每个点的位置，然后开始下一个周期的计算。
|
1 |
for (int i = 0 ; i < nBodies; i++) { //根据速度，计算新位置 p[i].x += p[i].vx*dt; p[i].y += p[i].vy*dt; p[i].z += p[i].vz*dt; } |
CUDA并行程序设计
并行的基本思路
由于在计算每一个点的受力情况是相互独立互不影响的，可以利用GPU高并发的特点，用n个线程并发计算速度变化，然后用n个线程并发更新位置。这样可以很大程度的提升代码的效率。
并行的详细设计
Step1：申请CPU和GPU内存空间并对数据进行初始化和拷贝操作。
|
1 |
//申请空间并初始化 int bytes = nBodies * sizeof(Body); float *buf; cudaMallocHost((void **)&buf,bytes); randomizeBodies(buf, 6 * nBodies); float *d_buf; cudaMalloc((void **)&d_buf,bytes); Body *d_p=(Body *)d_buf; cudaMemcpy(d_buf,buf,bytes,cudaMemcpyHostToDevice); |
Step2：设计bodyForce函数
对于每个线程，只需根据threadIdx.x和blockIdx.x计算对应在p中的位置，然后进行全局更新即可。
|
1 |
//并行的bodyForce函数 __global__ void bodyForce(Body *p, float dt, int n) { int i=threadIdx.x+blockIdx.x*blockDim.x; if(i<n){ float Fx = 0.0f; float Fy = 0.0f; float Fz = 0.0f; for (int j = 0; j < n; j++) { float dx = p[j].x - p[i].x; float dy = p[j].y - p[i].y; float dz = p[j].z - p[i].z; float distSqr = dx*dx + dy*dy + dz*dz + SOFTENING; float invDist = rsqrtf(distSqr); float invDist3 = invDist * invDist * invDist; Fx += dx * invDist3; Fy += dy * invDist3; Fz += dz * invDist3; } p[i].vx += dt*Fx; p[i].vy += dt*Fy; p[i].vz += dt*Fz; } } |
Step3：设计integrate_position函数
对于每个线程，只需根据threadIdx.x和blockIdx.x计算对应在p中的位置，然后更新位置即可。
|
1 |
__global__ void integrate_position(Body *p,float dt,int n){ int i=threadIdx.x+blockIdx.x*blockDim.x; if(i<n){ // integrate position p[i].x += p[i].vx*dt; p[i].y += p[i].vy*dt; p[i].z += p[i].vz*dt; } } |
优化思路
分析基本并行程序可知，程序的主要开销是在bodyForce中。对于每一个线程，需要计算的量是n，而且对于全局内存的频繁访问会导致大量的时间开销，故，我们可以用空间换取时间，使用BLOCK_STEP倍的线程使每个线程的计算量得到提升，同时在同一个线程块中使用shared_memory存放最近会访问到的数据。
空间优化
• 对于每个点，每次计算都会访问全局的数据，这会造成大量开销
• 对于一个线程块，块内有共享的内存，可以使用该内存每次获取一部分进行计算，减少访存 时间
空间换时间
• 引入 BLOCK_STEP 变量，每个线程仅计算一个点的部分数据
• 每个线程只执行循环 (n=4096) 被 BLOCK_STEP 划分后的一部分
下图展示了划分的原理，可以很好的理解优化的思路，下图中的n所在的每个块都有BLOCK_SIZE个数据，被BLOCK_STEP划分。注意理解数据划分和线程的联系。
优化2—— 计算合并
每个线程的计算时间不一致，在上述的程序逻辑中，需要等待所有线程计算完速度的改变量才进行位置的全局更新，我们可以将更新位置的操作放在计算速度改变量中，用计数器判断某个点的数据是否已经被算完，如果算完则立即更新，无需等待。Flag数组是用于记录计算情况的数据，每个线程算完后会对对应的flag[i]原子减一，如果为0则立即更新位置并重置flag[i]，注意操作的原子独立性。
|
1 |
atomicSub(&flag[i], 1); if(!atomicMax(&flag[i], 0)){ // integrate position atomicAdd(&p[i].x,p[i].vx*dt); atomicAdd(&p[i].y,p[i].vy*dt); atomicAdd(&p[i].z,p[i].vz*dt); atomicExch(&flag[i],BLOCK_STEP); } |
优化3—— 编译优化
使用编译优化，加入循环展开，调整得到最佳参数。
优化4—— 其他优化方向
考虑到GPU的硬件特点，一定是最后一组最后算完，故在最后一组顺便计算结果。
效果对比
串行效率
未经优化的并行
优化后的并行
效率有1771倍的提升，优化很成功。
其他思路&源码资源
部分代码改汇编也可以得到不错的提升。
还可以考虑用多个核来计算，这样会有差不多一个数量级的提升，当数据量增大时会更明显。
https://github.com/nightmare2423/cuda
https://github.com/nightmare2423/cuda
微信扫一扫打赏
支付宝扫一扫打赏
