大规模三维空间纤维距离计算程序[CUDA校园程序设计大赛]

很早就报名了CUDA校园程序设计大赛，一直就因为各种事情加上各种懒没动参赛作品的脑筋。今年的工作人员非常尽责啊，没提交作品打电话询问过两三次了，最近终于抽出完整的n个小时整理出了一个参赛作品，作品来源是为解决实验室研究过程中实际遇到的一个问题而写的cuda程序，并没有什么高级的并行改造之类的动作，就是充分利用SIMD特性而已哈。程序没有做什么改动，只是做了计时测试，写了个自认为还算认真的说明文档，倒是为了做对比而写的cpu程序遇到了一些有意思的事情。

作品的压缩包下载： http://www.opendrive.com/files/7150902_XFTG5_1f21/cuda.zip
以下便是说明文档的内容：
1. 问题描述
在医学可视化领域，纤维（可理解为空间曲线）经常被用于描述一些软组织器官结构，如肌肉、脑白质等。其中每个纤维都可由三维空间中的若干个控制点描述。而纤维之间的距离计算往往是各种分析处理工作的必须步骤。
若给定两个纤维集合FiberSetA与FiberSetB，其中FiberSetA中的纤维数为CntA，FiberSetB中的纤维数为CntB。我们的问题即为求解两纤维集合中两两纤维间的距离（即结果中共包含有CntA * CntB个距离）。
两空间曲线间的距离可以有许多种定义方式，本程序选用的是控制点到曲线间的距离求平均的方式。例如，对于纤维X（由x0, x1, …, xt-1控制点构成）和纤维Y（由y0, y1, …, ys-1控制点构成），设任意两个控制点x与y之间的距离为Dist(x,y)，则控制点x到纤维Y的距离定义为Dist(x,Y) = min{Dist(x,y0), Dist(x,y1), …, Dist(x,ys)}，相应的纤维X到纤维Y的距离定义为Dist(X,Y) =∑(xi, Y) / t。
当涉及的纤维数量和控制点较多时，传统的CPU计算程序所需时间会很长。

2. 算法与加速
该问题所对应的是计算密集型的任务（即求解大量的“距离”），因此在算法上并没有太大的优化空间。而在“增强并行性”上做文章无疑会带来显著的加速效果，因为不同距离之间的计算是完全不存在依赖性的。可以说，这是一个GPU拥有绝对优势的任务场合。
本程序中，每个thread（CUDA中的thread概念）对应一个距离。在只使用global memory的最初版本中就已经展现了较大的加速效果。
进一步探究可以发现，对于同一个纤维所涉及到的数据，每个与此纤维相关的thread都会从Global Memory中读取该纤维的相关数据，而访问Global Memory的延迟是比较高的，对于处于一个block中且涉及相同纤维的thread而言，我们可以利用访问延迟比Global Memory低很多的Shared Memory存储最初版本中重复读取的数据来实现进一步的加速。实践证明，这样的优化再次带来了可观的加速效果。（注：当然，在纤维的控制点太多时，由于Shared Memory 16kByte的限制，并不能将所有最初版本中重复读取的数据都放进来，这也是目前本程序的局限性之一。）

3. 编程与运行环境
本程序采用VS2008(MSVC9) + CUDA Toolkit 3.1 + CUDA SDK 3.1开发；
对比的CPU程序分别通过VS2008(MSVC9)、VS2010(MSVC10)和Mingw(GCC4.5.1)编译，其中MSVC9和MSVC10的主要编译选项包括：/O2、/MT；而GCC4.5.1的主要编译选项包括：-s、-O3。皆编译为x64平台程序。同时参赛提交的压缩包中“bin/x86”文件夹中也放入了适用于32位平台的程序。（注：CPU程序并没有做如使用指令集、多线程等的特意优化。）
测试运行环境为：Windows7(64bit) + Intel E7500 + 4G DDR2 800 + GTX260

4. 测试结果
重点终于来了。每个程序运行3次取计时结果的平均值，单位为秒，精确到小数点后第三位，分别针对数据类型为float和double进行了测试。测试数据为data文件夹下“a.txt”和“b.txt”描述的两个纤维集，其中a和b中分别含有5265个和360个纤维，控制点个数都为312个。测试结果统计见下表：

	CUDA3.1	MSVC9	MSVC10	GCC4.5.1
float	1.173s	242.582s	64.834s	55.799s
double	5.510s	237.620s	65.306s	55.523s

下图更加直观的体现了速度上的差异（没有将过慢的MSVC9编译版本列入）：



有惊喜，也有惊讶。惊喜的是针对本任务的计算场合，CUDA的确可以带来相当可观的加速效果；惊讶的是，同是CPU版本，对于同样的代码MSVC9和GCC4.5.1编译出的程序表现差距居然如此之大，而且还是在Microsoft自家的Windows平台上，真是令人匪夷所思（我已仔细检查过VS2008的工程设置，确定没有使用Debug编译J）。扯远了，欢迎有兴趣的朋友邮件交流此事。

5. 总结与未来工作
本程序并不涉及“并行改造”这样往往很复杂的算法优化工作，因为其本身面临的就是一个“可并行计算密集型”的任务，与CUDA架构的优势之处不谋而合，产生可观的加速效果也就在情理之中了。
当然，本程序在开发的时候并没有过多考虑通用性和扩展性，目前还存在以下几点问题：
1) 前文已经提到的Shared Memory容量限制的问题。本程序默认Shared Memory不会发生不够用的情况，并没有考虑超出容量限制时的处理措施；
2) 对比的CPU程序还存在很大的优化空间：如使用SSE指令集，使用多线程等；
3) 当Kernel函数执行超过一定时间时（Win7/Vista下3秒，WinXP下6秒），Windows显示设备检测机制会中断GPU任务的执行。本程序只采用简单地将原Grid等分为固定的N份然后多次执行Kernel的方法，其中N为程序中人为定义的常数。若根据数据规模和当前运行环境的设备性能参数自动调节这个N会使程序更加具备通用性。

6. Bonus(Tips)
1) 假设一个a.h文件被一个b.cu包含，在b.cu被build之后，如果a.h的内容发生改动再去build b.cu，nvcc并不会发觉a.h的内容已经改变从而跳过对b.cu文件的编译，这个可大可小的bug如果忽视很容易带来很多无谓的调试工作；
2) 关于上文提到的Kernel函数执行时间过长所遇到的问题，有一个治标不治本的方法，就是修改注册表，详情请见：
http://www.microsoft.com/whdc/device/display/wddm_timeout.mspx
不过该方法并不推荐，权且当作应急之需吧。
3) 关于VS中的Code Definition Window：个人非常喜欢用VS中的Code Definition Window来看代码，不过它只支持VS认识的文件格式（.c/.cpp/.h等等），由于VS并不认识.cu文件，因此Code Definition Window就不能用来快速了解.cu文件中的变量和函数信息了。一个解决方法是，新建一个temp.cpp文件，并将其设置为编译范围之外，然后将.cu文件的内容拷贝过去，这样在浏览这个temp.cpp文件的时候Code Definition Window就又开始工作了。