多核CPU加速并行计算的快捷开发和应用

　　不知有多久没来CSDN写文章了，今天终于勤快了一回，把最近一周的一点业余研究心得总结成这个文章。研究这个问题，是为了寻找在OpenCL之外，可以优化CPU执行效率的方法，优化对三维场景中大量物体实时计算的效率。比如拾取，变换，物理效果之类。之前实践了OpenCL和DirectCompute，但是只能在公司的高级电脑上，在很多只有集成显卡而且还跑着WinXP的电脑上是无缘这些了。而本文要介绍的这个方法实现起来比较简单，只要有多核CPU就能加速，对操作系统也几乎没要求，更不需要搞复杂的
CUDA, OpenCL，DirectCompute 之类，所以我这种又懒又穷的人就自己用了。

设计思想

　　利用多核CPU同时执行多个线程，使每个线程运行在一个CPU核心上，让线程保持连续多次运行，实现了并行计算。整个功能只用到了 Windows API，实现代码极少，根据这种实现思路，只要在提供线程API的操作系统中都可以很容易的实现。对于N核CPU可以达到N倍的并行加速效果。当然，这个用法只限于使用线程函数，加速效果也不能和显卡相比。

实现方法

　　功能的实现主要是两部分，线程池和线程并行控制。线程池提高了线程连续运行和重复使用的效率，线程并行控制最终实现了并行加速。

　　(1).线程池 (Thread Pool)

　　线程池是用来维持线程生命周期，重复执行多次任务，减少操作系统维护时间的一种线程管理方法。一个线程在创建后，一旦运行结束，这个线程就不存在了，操作系统在创建线程和线程执行结束时都需要做很多底层操作。如果执行大量的异步任务，只是简单的给每次任务逐个创建线程，操作系统会浪费很多时间和资源，最终多个线程并行的时间很难比顺序执行更快。线程池使一个线程重复执行多次线程函数，而这个线程不会立刻结束，这样就保证了线程在运行多个任务时始终存在。

　　最简单的线程池的实现如下：

//用户函数类型：
typedef void (*UserThreadProc)();
//创建一个执行队列(其实就是数组)，把它作为参数传递给下面的内部线程函数，在线程中执行这个队列。
std::vector<UserThreadProc> tasklist;
//内部的线程函数，这个是线程直接运行的函数，在这个函数里调用队列中的函数。
DWORD __stdcall  _InternalThreadProc(LPVOID args)
{
std::vector<UserThreadProc>* tasklist = (std::vector<UserThreadProc>*)args;
for(std::vector::iterator i = tasklist->begin(); i!=tasklist->end(); i++)
{
(*i)();
}
tasklist->clear();
return NULL;
};
//应用程序如下:
void MyTaskFunc()//实现一个线程函数作为任务。
{
//Do some hard work.
};
void main()
{
// 向队列中放入100个任务，这里只是简单的重复调用”MyTaskFunc()”。
tasklist.resize(100, MyTaskFunc);
// 使用Windows API创建一个线程，但是指定线程函数为"_InternalThreadProc"，把执行队列"tasklist"作为参数传入。
HANDLE thread = CreateThread(NULL, NULL, _InternalThreadProc, &tasklist, NULL, NULL);
// 现在这个线程开始执行这100次任务了。
WaitForSingleObject(thread, INFINITE);//等待完成即可。
};

　　这个线程池只是个简单的例子，没有用户参数，也不能控制线程的运行，毫无实用意义，需要继续改进。改进后如下：

//定义一个有参数和返回值的用户函数类型
typedef __int64 (*UserThreadProc)(__int64 user_args);
//然后定义一个结构保存异步执行的信息。
struct AsyncState
{
UserThreadProc func; //用户执行的函数。
__int64 params; //用户函数的参数。通过保存地址可以转换成任意类型。
__int64 result; //用户函数的返回值。通过保存地址可以转换成任意类型。
bool is_finished; //指示是否执行完成。由底层来设置，调用者就别随意修改了。
};
//现在可以定义一个线程池对象：
class CThreadPool
{
private:
bool                        m_enable;
HANDLE                      m_thread;
std::queue<AsyncState*>     m_tasklist;
public:
CThreadPool();
~CThreadPool();
void enqueue(AsyncState* user_call);
size_t curr_pending_tasks();
};
//然后实现一个改进的线程函数，把线程池对象自身当做参数传入。
DWORD __stdcall _InternalThreadProc(LPVOID args)
{
CThreadPool* tp = (CThreadPool*)args;
while(tp->m_enable)
{
while(tp->m_tasklist.size())
{
AsyncState* as = tp->m_tasklist.front();
as->result = as->func(as->params);
as->is_finished = true;
tp->m_tasklist.pop();
}
SuspendThread(tp->m_thread);
}
};
//然后CThreadPool可以用以下的代码实现：
CThreadPool::CThreadPool()
{
this->m_enable = true;
this->m_thread = CreateThread(NULL,NULL,_InternalThreadProc,this,CREATE_SUSPENDED,NULL);
}
CThreadPool::~CThreadPool()
{
this->m_enable=false;
ResumeThread(this->m_thread);
WaitForSingleObject(this->m_thread,INFINITE);
}
void CThreadPool::enqueue(AsyncState* user_call)
{
user_call->is_finished=false;
this->m_tasklist.push(user_call);
if (this->m_tasklist.size()==1)
{
ResumeThread(this->m_thread);
}
}
size_t CThreadPool::curr_pending_tasks()
{
return this->m_tasklist.size();
}

　　把_InternalThreadProc()和CThreadPool的实现代码放在cpp文件中，这样_InternalThreadProc()作为内部使用的函数，使用者是看不到的，使用者只需设置好所有的AsyncState，通过CThreadPool::enqueue()加入队列，这些函数就可以在线程中自动运行了。创建线程时使用了"CREATE_SUSPENDED"，线程不会立刻运行，每当检测到初次加入了AsyncState，通过"ResumeThread(this->m_thread);"唤醒线程，当所有用户函数都执行完毕时，线程被挂起，"SuspendThread(tp->m_thread);"，直到被再次唤醒。只有析构函数中的代码才能让线程结束。

　　(2).线程并行控制

　　线程并行控制的目的是，当线程池中包含了多个线程时，均衡多个线程和CPU核心的任务，让每个线程运行在最佳状态，相互影响最小，总的运行时间最短。这部分是实现并行的重点。具体功能是两个，一个是分配线程核心，一个是均分线程任务。

　　分配线程核心：利用多核CPU的天然并行能力，使每个线程运行运行在一个核心上，每个线程独享一个核心的计算能力。Windows系统上的实现方法很简单，直接使用"SetThreadAffinityMask()"，每个核心独立运行互不影响时，执行效率比较高。

　　均分线程任务：在N核心CPU的系统中，线程池可以有N个线程(这是最佳用法)，添加异步任务时，应该始终把任务加入到等待队列最小的一个线程(通过检查"CThreadPool::curr_pending_tasks()")，使完成任务的总时间最少，这个简单的分配方法有效的减少了线程同步花费的时间。

测试代码

　　根据上述内容，实现了一个比较通用的并行调用接口，封装了以上所有功能，线程数量自动和CPU核心数量一致。在接口层面，调用者只需设置多个AsyncState，依次通过enqueue传入即可。

　　测试内容是生成10,000,000个随机位置的 Axis Aligned Bounding Ellipse 并进行射线求交测试，直接遍历计算，不使用八叉树等空间优化方法。并行计算时，程序根据CPU核心数自动创建4个线程，把10,000,000次测试分成4路同时执行，每个线程执行2,500,000次计算，并在计时中加入了线程同步的时间。说明一下，这只是一个夸张的测试例子，实际场景中不可能有这么多对象，而且场景管理肯定需要空间八叉树来优化。

　　这里附上射线与椭球求交算法的主要代码：

bool Ray3D::IsIntersectEllipsoid(Region3D& obj, float* result_distance, Vector3* result_position) const
{
//计算方法：

//椭球标准方程：
// ((x - cx) / rx)^2 + ((y - cy) / ry)^2 + ((z - cz) / rz)^2 = 1
//射线参数方程：
// x = bx + dx * t；y = by + dy * t；z = bz + dz * t；
//产生关于参数t的交点方程：
// ((dx * t + bx - cx) / rx)^2 + ((dy * t + by - cy) / ry)^2 + ((dz * t + bz - z) / rz)^2 = 1
//化简形式为：
// a * (t^2) + b * t + c = 0
// a = (Ry*Rz*Dx)^2+(Rx*Rz*Dy)^2+(Rx*Ry*Dz)^2
// b = ((Ry*Rz)^2*Dx*(Bx-Cx)+(Rx*Rz)^2*Dy*(By-Cy)+(Rx*Ry)^2*Dz*(Bz-Cz))*2
// c = (Ry*Rz*(Bx-Cx))^2+(Rx*Rz*(By-Cy))^2+(Rx*Ry*(Bz-Cz))^2-(Rx*Ry*Rz)^2
//算出t即可。

//实现过程：
float RxRy2 = obj.SizeX() * obj.SizeX() * obj.SizeY() * obj.SizeY() * 0.0625f;
float RxRz2 = obj.SizeX() * obj.SizeX() * obj.SizeZ() * obj.SizeZ() * 0.0625f;
float RyRz2 = obj.SizeY() * obj.SizeY() * obj.SizeZ() * obj.SizeZ() * 0.0625f;
float RxRyRz = obj.SizeX() * obj.SizeY() * obj.SizeZ() * 0.125f;
float Bx_Cx = this->PosX - obj.Center().x;
float By_Cy = this->PosY - obj.Center().y;
float Bz_Cz = this->PosZ - obj.Center().z;
float a = RyRz2 * this->DirX * this->DirX + RxRz2 * this->DirY * this->DirY + RxRy2 * this->DirZ * this->DirZ;
float b = (RyRz2 * this->DirX * Bx_Cx + RxRz2 * this->DirY * By_Cy + RxRy2 * this->DirZ * Bz_Cz) * 2.0f;
float c = (RyRz2 * Bx_Cx * Bx_Cx + RxRz2 * By_Cy * By_Cy + RxRy2 * Bz_Cz * Bz_Cz) - RxRyRz * RxRyRz;
float d = b * b - a * c * 4.0f;
if(d < 0.0f)
return false;
float r_distance = (-b - sqrtf(d)) / a * 0.5f;
Vector3 r_position = Vector3(
this->PosX + this->DirX * r_distance,
this->PosY + this->DirY * r_distance,
this->PosZ + this->DirZ * r_distance);
if(result_distance != 0)
(*result_distance) = r_distance;
if(result_position != 0)
(*result_position) = r_position;
return true;
}

测试环境1: Intel Xeon E3 1225, 3.1 GHz, 4 Cores (No HT), Windows 7 x64.

以下是运行5次的结果：

1:


2:



3:



4:



5:



测试环境2: Intel Core2 Duo T5800, 2.0 GHz,2 Cores (No HT), Windows XP x64.

以下是运行5次的结果：

1:



2:



3:



4:



5:



　　结果还算理想。
　　不过在第一个测试中，并行运行的加速倍率经常可以超过N核心数(4倍以上)，不知计时方法是否有问题。Xeon E31225 没有睿频加速，可以排除这个影响，计时用的是QueryPerformanceFrequency()和QueryPerformanceCounter()，也够精确了啊。还求了解CPU的大牛指点一下。