Cuda Nbody笔记

Nbody Cuda ： N个元素，切割成N/p个块，每个块p个元素。
开设N/p个block,每个block里面开p个线程。
每个线程计算一个元素和其他N个元素的受力和加速度。
每个线程将需要计算的N个元素分成N/p次完成，每次计算p个元素。
calculate_forces：总的入口，计算一个元素和所有元素的受力和加速度。
tile_calculation：线程一次计算和p个元素的受力。
bodyBodyInteraction：两个元素之间的受力和加速度。

程序中的tile是指p*p个元素的方块。

受力公式如下：万有引力公式

 

 

__device__ float3 bodyBodyInteraction(float4 bi, float4 bj, float3 ai) { float3 r; // r_ij [3 FLOPS] r.x = bj.x - bi.x; r.y = bj.y - bi.y; r.z = bj.z - bi.z; // distSqr = dot(r_ij, r_ij) + EPS^2 [6 FLOPS] float distSqr = r.x * r.x + r.y * r.y + r.z * r.z + EPS2; // invDistCube =1/distSqr^(3/2) [4 FLOPS (2 mul, 1 sqrt, 1 inv)] float distSixth = distSqr * distSqr * distSqr; float invDistCube = 1.0f/sqrtf(distSixth); // s = m_j * invDistCube [1 FLOP] float s = bj.w * invDistCube; // a_i = a_i + s * r_ij [6 FLOPS] ai.x += r.x * s; ai.y += r.y * s; ai.z += r.z * s; return ai; }

__device__ float3 tile_calculation(float4 myPosition, float3 accel) { int i; extern __shared__ float4[] shPosition; //each thread will process blockDim.x elements. for (i = 0; i < blockDim.x; i++) { accel = bodyBodyInteraction(myPosition, shPosition[i], accel); } return accel; }

__global__ void calculate_forces(void *devX, void *devA) { extern __shared__ float4[] shPosition; float4 *globalX = (float4 *)devX; float4 *globalA = (float4 *)devA; float4 myPosition; int i, tile; float3 acc = {0.0f, 0.0f, 0.0f}; int gtid = blockIdx.x * blockDim.x + threadIdx.x; //current element to be processed myPosition = globalX[gtid]; //N elements. p elements are processed parallelly for (i = 0, tile = 0; i < N; i += p, tile++) { //get element in the tile int idx = tile * blockDim.x + threadIdx.x; shPosition[threadIdx.x] = globalX[idx]; //synchronize to make sure that all the elements in the same tile has been fetched into shared memory. __syncthreads(); acc = tile_calculation(myPosition, acc); //synchronize to make sure all the threads have been finished calculation. __syncthreads(); } // Save the result in global memory for theintegration step. float4 acc4 = {acc.x, acc.y, acc.z, 0.0f}; globalA[gtid] = acc4;//save the acceleration }