14.7 CPU的优化

在CUDA移植的论文中，会经常和没有进行最优化处理的CPU实现相比较。虽然在这些论文中描述的负载上，CUDA硬件的计算速度比CPU快，但是得到的加速比往往比一旦对CPU实现予以优化的结果高。

为了得到一些有助于在GPU与最先进的CPU优化间进行权衡的启示，我们将使用让多核CPU达到峰值性能必要的两个关键策略来优化N-体计算。

·SIMD指令可以在一个指令周期执行多个单精度浮点数的计算。

·当CPU有可用的执行核心时，多线程将实现趋近于线性的加速比。从2006年开始，多核CPU已经被广泛地使用，预计N-体计算的加速比和核心数量呈线性关系。

既然N-体计算有如此高的计算密度，我们不必关注亲和性（例如，尝试使用NUMA API来关联内存缓冲区到特定CPU）。这个计算存在大量重用数据，CPU的缓存将使得计算仅需少量内存的流量。

20世纪90年代末期，从奔腾III（PentiumIII）开始，英特尔的x86构架增加了SIMD流扩展指令集（SSE）。它们增加了一组八个128位的

XMM寄存器，能够操作4个被捆绑在一起的32位浮点数。[1]例如，ADDPS指令使用XMM寄存器能实现4个浮点数的并行相加。

当将N-体计算采用SSE指令集，我们一直在使用的AOS（结构体数组）内存结构将成为一个问题。虽然个体数据和XMM寄存器位数一样是16个字节，指令集却要求我们重新排列数据，将x、y、z和mass分量放在独立的寄存器。而不是在计算个体间作用力时执行这一操作，我们重新排列内存结构为“数组结构”：不是使用一个float4数组（每个元素包含给定个体的x、y、z和mass），我们使用4个float数组、即x坐标数组、y坐标数组、z坐标数组和mass坐标数组。数据重新排列后，仅通过4个机器指令即可将4个个体的数据载入XMM寄存器；4个个体位置的差向量将通过3条SUBPS指令计算得到；等等。

为了简化SSE指令集代码编写工作，英特尔与编译器厂商合作，使SSE指令集获得跨平台的支持。特殊的数据类型_m128对应于128位的寄存器和操作数大小，可以用于内置函数，例如对应于SUBPS指令的_mm_sub_ps()。

为了便于N-体计算的实现，我们还需要使用全精度的平方根倒数实现。SSE指令集有一个RSRTPS内置指令来计算近似的平方根倒数，但是它的12位的估计值必须使用牛顿-拉夫逊迭代精化到完整浮点精度。[2]

代码清单14-8给出了SSE指令集在计算个体间作用力的实现，它带有2个_m128类型变量代表2个个体数据，并行地计算4个个体间的作用力，并返回3个最终力向量。代码清单14-8描述的功能等同于代码清单14-1和14-2，但可读性明显减弱。注意：x0、y0和z0变量包含同一个个体的数据，通过复制4次_m128型变量得到。

为了利用多核的优势，我们必须使用多个线程，并且让每一个线程承担部分计算任务。多核CPU使用和多GPU相同的策略：[3] 只需要把输出行平均划分给每一个线程（每个CPU核一个线程），在每一时间步长，父线程通过信号指挥工作线程进行工作并等待工作线程结束。因为线程的创建需要时间并且可能失败，我们的应用在初始化时创建一个CPU线程池，并使用线程同步让工作线程等待它们的工作，当工作完成时发出信号。

本书的可移植线程库，详见附录A-2。库中实现了函数 processorCount()，返回系统中CPU核的数目，并有一个C++类 wakerThread包含了创建和销毁CPU线程以及以同步或异步的方法委派工作的功能。在使用delegateAsynchronous()成员函数异步委派工作后，调用静态函数waitAll()等待工作线程的结束。

代码清单14-9给出了分配N-体计算任务到CPU工作线程的代码。sseDelegation结构是用来向每一个工作线程传达委派信息的；delegateSynchronous函数调用需要传递一个函数指针和一个将传递给传入函数的void类型的指针（在这个例子中，void指向对应CPU线程的sseDelegation结构）。

代码清单14-8 个体间相互作用（SSE版本）

static inline _m128   
rcp_sqrt(nr_ps(const _m128 x)   
{ const _m128 nr  $=$  _mm_rsqrt_ps(x), muls  $=$  _mm_mul_ps(_mm_mul_ps(nr, nr), x), beta  $=$  _mm_mul_ps(_mm_set_psl(0.5f), nr), gamma  $=$  _mm_sub_ps(_mm_set_psl(3.0f), muls); return _mm_mul_ps(beta, gamma);   
}   
static inline _m128 horizontal_sum_ps(const _m128 x)   
{ const _m128 t  $=$  _mm_add_ps(x, _mm_movehl_ps(x, x)); return _mm_add_ss(t, _mm_shuffle_ps(t, t, 1));   
}   
inline void bodyBodyInteraction( _m128& f0, _m128& f1, _m128& f2, const _m128& x0, const _m128& y0, const _m128& z0, const _m128& x1, const _m128& y1, const _m128& z1, const _m128& massl,

const _ml28& softeningSquared)   
{ _ml28 dx  $=$  _mm_sub_ps(x1, x0); _ml28 dy  $=$  _mm_sub_ps(y1, y0); _ml28 dz  $=$  _mm_sub_ps(z1, z0); _ml28 distSq = _mm_add_ps( _mm_add_ps( _mm.mul_ps(dx, dx)), _mm.mul_ps(dy, dy)) , _mm.mul_ps(dz, dz) ); distSq  $=$  _mm_add_ps(distSq, softeningSquared); _ml28 invDist  $=$  rcp_sqrt(nr_ps(distSq)); _ml28 invDistCube  $=$  _mm.mul_ps( invDist, _mm.mul_ps( invDist, invDist ) ）; _ml28 s  $=$  _mm.mul_ps(mass1, invDistCube); f0  $=$  _mm_add_ps(a0, _mm.mul_ps(dx, s)); f1  $=$  _mm_add_ps(al, _mm.mul_ps(dy, s)); f2  $=$  _mm_add_ps(a2, _mm.mul_ps(dz, s));

代码清单14-9 多线程SSE指令集（主控线程代码）

float   
ComputeGravitation_SSE_threaded( float \*force[3], float \*pos[4], float \*mass, float softeningSquared, size_t N   
} chTimerTimestamp start, end; chTimerGetTime( &start ); { sseDelegation \*psse  $=$  new sseDelegation[g_numCPUCores]; size_t bodiesPerCore  $= \mathrm{N}$  /g_numCPUCores; if(N  $\text{十}$  g_numCPUCores）{ return 0.of; } for（size_t i = 0;i < g_numCPUCores;  $\mathbf{i + + }$  ）{ psse[i].hostPosSOA[0]  $=$  pos[0]; psse[i].hostPosSOA[1]  $=$  pos[1]; psse[i].hostPosSOA[2]  $=$  pos[2]; psse[i].hostMassSOA  $=$  mass; psse[i].hostForceSOA[0]  $=$  force[0]; psse[i].hostForceSOA[1]  $=$  force[1]; psse[i].hostForceSOA[2]  $=$  force[2]; psse[i].softeningSquared  $=$  softeningSquared; psse[i].i  $=$  bodiesPerCore\*i; psse[i].n  $=$  bodiesPerCore; psse[i].N  $=$  N;

g_CPUThreadPool[i].delegateAsynchronous( sseWorkerThread, &psse[i]);   
} workerThread::waitAll(g_CPUThreadPool,g_numCPUcores）; delete[] psse;   
} chTimerGetTime（&end）; return(float)chTimerElapsedTime(&start,&end）\*1000.of;

最后，代码清单14-10给出了sseDelegation结构和被代码清单14-9中ComputeGravitation_SSE_threaded调用的委派函数。它实现一次完成4个个体间作用力的计算，在存储最后的输出结果前使用horizontal_sum_ps()函数累计得到的4个部分和。这个函数以及其调用的所有函数，使用SOA的内存结构进行所有的输入和输出。

代码清单14-10 从属线程代码sseWokerThread

struct sseDelegation {
    size_t i; // base offset for this thread to process
    size_t n; // size of this thread's problem
    size_t N; // total number of bodies
    float *hostPosSOA[3];
    float *hostMassSOA;
};

float *hostForceSOA[3]; float softeningSquared;   
}；   
void sseWorkerThread( void \*p) { sseDelegation \*p = (sseDelegation \*) p; for (int k = 0; k < p->n; k++) { int i = p->i + k; _m128 ax = _mm_setzero_ps(); _m128 ay = _mm_setzero_ps(); _m128 az = _mm_setzero_ps(); _m128 \*px = (_m128 \*) p->hostPosSOA[0]; _m128 \*py = (_m128 \*) p->hostPosSOA[1]; _m128 \*pz = (_m128 \*) p->hostPosSOA[2]; _m128 \*pmass = (_m128 \*) p->hostMassSOA; _m128 x0 = _mm_set_psl{ p->hostPosSOA[0][i] }； _m128 y0 = _mm_set_psl{ p->hostPosSOA[1][i] }； _m128 z0 = _mm_set_psl{ p->hostPosSOA[2][i] }； for (int j = 0; j < p->N/4; j++) { bodyBodyInteraction{ ax, ay, az, x0, y0, z0, px[j], py[j], pz[j], pmass[j], _mm_set_psl{ p->softeningSquared }); } // Accumulate sum of four floats in the SSE register ax = horizontal_sum_ps(ax); ay = horizontal_sum_ps(ay); az = horizontal_sum_ps(az); _mmstore_ss( float \*) &p->hostForceSOA[0][i], ax); _mmstore_ss( float \*) &p->hostForceSOA[1][i], ay); _mmstore_ss( float \*) &p->hostForceSOA[2][i], az); }

[1] 英特尔后来又增加了一些指令，能够让XMM寄存器捆绑处理多个整型数（多达16字节）或者两个双精度的浮点数，但是我们没有使用这些特性。我们也没有使用AVX（高级向量扩展）指令集。AVX的特性在于寄存器和指令支持两倍位宽的单指令多数据操作（256位），因此它可能拥有两倍于当前性能的潜能。
[2] 这段代码没有出现在SSE编译器支持中而且特别难找。我们的实现来自http://nume.google.com/svn/trunk/fosh/src/sse_approx.h。

[3] 事实上，在第9章的多线程多GPU中使用了相同的独立于平台的线程库。