14.2 简单实现

代码清单14-1给出了上一节描述的个体之间相互作用的函数实现，使用了__host__和__device__关键词双重修饰，CUDA编译器就会知道它既可以被CPU也可以被GPU调用。这个函数带有模板，因此它可以被float和double同时使用（在这本书中，只有float型是完全实现了的）。它返回三维作用力矢量，用三元组（ax,ay,az）[1]表示。

代码清单14-1 bodyBodyInteraction函数

template <typename T>
    __host__device__void bodyBodyInteraction(
        T& ax, T& ay, T& az,
        T x0, T y0, T z0,
        T x1, T y1, T z1, T mass1,
        T softeningSquared)
{
    T dx = x1 - x0;
    T dy = y1 - y0;
    T dz = z1 - z0;
    T distSqr = dx * dx + dy * dy + dz * dz;
    distSqr += softeningSquared;
    T invDist = (T)1.0 / (T)sqrt(distSqr);
    T invDistCube = invDist * invDist * invDist;
    T s = mass1 * invDistCube;
    ax = dx * s;
    ay = dy * s;
    az = dz * s;
}

代码清单14-2给出了计算每个个体受到所有作用力的函数实现。对于每个个体，它把自己的位置加载到（myX, myY, myZ），然后调用函数bodyBodyInteraction来计算和其他每个个体之间的作用

力。函数名中的AOS，表示输入的数据来自一个结构体数组（array of structures，AOS）。包含四个float型数的结构体（x,y,z,mass）指定了个体的位置和质量。使用float4类型的数据是为了方便GPU的实现，这样可以得到硬件在加载和存储上的天然支持。而我们优化的CPU实现，详见第14.9小节，使用包含x、y、z和mass四个float型数组组成的结构数组（structure of arrays，SOA），这种形式能够有利于SIMD（单指令多数据）指令集的处理。结构数组不是很适合GPU实现，因为其中所需的四个基指针会耗费过多的寄存器。

代码清单14-2 ComputerGravitation_AOS函数（CPU实现）

float   
ComputeGravitation_AOS( float \*force, float \*posMass, float softeningSquared, size_t N   
1 chTimerTimestamp start, end; chTimerGetTime( &start ); for ( size_t i = 0; i < N; i++) { float ax  $=$  0.0f; float ay  $=$  0.0f; float az  $=$  0.0f; float myX  $\equiv$  posMass[i\*4+0]; float myY  $\equiv$  posMass[i\*4+1]; float myZ  $\equiv$  posMass[i\*4+2]; for ( size_t j = 0; j < N; j++) { float acc[3]; float bodyX  $\equiv$  posMass[j\*4+0];

float bodyY = posMass[j*4+1];  
float bodyZ = posMass[j*4+2];  
float bodyMass = posMass[j*4+3];  
bodyBodyInteraction<float>(ax, ay, az, myX, myY, myZ, bodyX, bodyY, bodyZ, bodyMass, softeningSquared);  
ax += acc[0];  
ay += acc[1];  
az += acc[2];  
}  
force[3*i+0] = ax;  
force[3*i+1] = ay;  
force[3*i+2] = az;  
}  
chTimerGetTime(&end);  
return (float) chTimerElapsedTime(&start, &end) * 1000.0f;

代码清单14-3给出了等价于代码清单14-2中函数的GPU版本。对于每一个个体，累计来自其他每一个个体的作用力值，然后将这个累计值写到一个作用力数组中。SM 2.x及以后版本中有一级和二级缓存，可以很好地加速这类负载，因为在内层循环中有大量的数据重用发生。

外层循环和内存循环都会将输入数组变成float4形式，这样可以保证编译器能够正确的发出单条16位加载指令。循环展开是在GPU上进行N-体计算最常用的优化方案，不难想象为什么是这样：GPU上分支结构的开销远大于CPU，因此在每次循环中减少指令能够带来很大好处，而且循环展开能够得到更多指令级并行（instruction level parallelism，ILP）的机会，这样GPU可以隐藏指令执行延时和内存访问延时。

表14-1 简单内核上的循环展开

为了在我们的N-体计算中利用循环展开的优点，我们需要在内层针对j的for循环前插入这行代码：

pragma unroll

不幸的是最佳循环展开因子需要依靠经验决定。表14-1总结了在这个内核中循环展开的效果。

在这个内核案例中，不使用循环展开时，每秒仅能计算个体间的作用力数目为250亿。即使展开因子是2，也能将计算性能提升到了300亿个每秒；将展开因子增加到16时，观测到内核的最佳性能：每秒计算个体间的作用力高达343亿个，得到 $37\%$ 的性能提升。

代码清单14-3 ComputeNBodyGravitation_GPU_AOS

template<typename T> global __void
ComputeNBodyGravitation_GPU_AOS()
{
    T *force,
    T *posMass,
    size_t N,
    T softeningSquared()
{
        for (int i = BlockIdx.x*blockDim.x + threadIdx.x;
            i < N;
            i += blockDim.x*gridDim.x)
        {
            T acc[3] = {0};
            float4 me = ((float4 *) posMass)[i];
            T myX = me.x;
            T myY = me.y;
            T myZ = me.z;
        for (int j = 0; j < N; j++) {
            float4 body = ((float4 *) posMass)[j];
            float fx, fy, fz;
            bodyBodyInteraction(
                &fx, &fy, &fz,
                myX, myY, myZ,
                body.x, body.y, body.z, body.w,
                softeningSquared);
            acc[0] += fx;
            acc[1] += fy;
            acc[2] += fz;
        }
        force[3*i+0] = acc[0];
        force[3*i+1] = acc[1];
        force[3*i+2] = acc[2];
}

[1] 原文误作（fx, fy, fz）。——译者注

14.2_简单实现

14.2 简单实现

代码清单14-2 ComputerGravitation_AOS函数（CPU实现）

表14-1 简单内核上的循环展开