优化无极限：从96秒到1.5秒
最近我们部门组织了一次编程竞赛。题目是这样的：在三维空间中，给出一组射线和一组三角形，其中每条射线给出基点O和方向量D，每个三角形给出三个顶点A,B,C的坐标，要求输出这组射线和这组三角形的所有交点。出题者会给出数据量很大的输入文件，可能包括几万条射线和几万个三角形，参加者提交自己的exe文件，在计算结果正确的前提下，运行时间最少者获胜。提交期限是一个星期。
其实这是光线跟踪研究领域的基本问题，基本算法已经比较成熟，所以上手比较简单。我和项目组的几个同事都跃跃欲试。说实话，我虽然做开发工作两年了，还确实没怎么研究过算法和优化。在为这个看似简单的算法题忙了5天后，感慨良多。我决定将这些感悟作为我在博客园开博后的首篇文章。
第一天：老牛拉破车的初始版本
琐事用百度，学术靠google。拿到这个问题后, 不用说，先 google一把，keyword就用“ray triangle intersections”,果然有很多相关论文。简单浏览了几篇文章，大致思想是<1>由三角形三个顶点得到三角形所在平面方程<2>求出射线和三角形所在平面的交点<3>判断该交点是否在三角形内部。射线可以用参数方程P(t)=O+tD表示，三角形平面方程可表示为Ax+By+Cz=0，联立方程可解出交点P。判断P是否在三角形ABC内，我用了一个比较简单的面积判别方法：分别求出面积S(PAB),S(PBC),S(PAC),三角形面积是S(ABC),如果S(PAB)+S(PBC)+S(PAC)<S(ABC),则P点在三角形内。花了两个小时写完后，用了一组样本测了一下，不错，结果正确。
第二天： 算法的精髓还是在数学
早上刚进公司，大家都在讨论这次的算法竞赛，似乎听说某人写的程序计算数据量最大那组文件（22600个三角形，18795条射线）只用了12秒。于是我兴致勃勃的用我的程序测了一下最大文件，受打击了，结果是96秒，惨不忍睹。与几个参加竞赛的同事讨论了一下，好像比较快的算法是用三角形重心坐标。什么是重心坐标？这个算法和三角形重心有什么关系？我只隐隐约约的记得在中学时学过三角形的重心是三角形三条垂线的交点。显然没有这么简单，发现自己理论知识储备不够啊。赶紧上Wiki百科，查了一下三角形重心坐标系，原来是这样的：如果空间一点P满足（u+v+w）P=uA+vB+wC,则u,v,w是点P在三角形ABC中的重心坐标。如果u+v+w=1（0<u,v,w<1）, 则点P在三角形内部。利用重心坐标，算法可以简化为用克拉默法则解行列式方程O+tD=(1-u-v)A+uB+vC，然后判断是否满足t>0,0<u,v<1。从算法上上看，这种方法比我第一天所用的方法计算量上明显小了很多。于是我按照此法写了第二个版本，果然同样那组文件运行时间降到48秒。虽然离最快的12秒还有距离，但较之第一个版本我已提高了很多，还是略有欣喜。
第三天：优化！优化！再优化！
再次认真的搜索查阅相关的资料后，我确信从理论上来说 ，几乎没有比这个还快的算法了。后面的工作只有两个字：优化！ 我一遍一遍的看着代码，哪里可以优化呢？我首先想到将频繁调用的函数换成宏。因为我知道任何函数调用都有入栈和出栈的开销，而宏会把代码直接插入到调用它的函数中，过多的用宏会增大可执行文件的尺寸，但由于没有出入栈的开销，运行速度会相对快些。但是快多少，我没有底。于是我把几个频繁调用的几个函数特别是在两重for循环里的一个函数换成宏，重新编译运行，提速效果然明显，又降了10秒多。优化了半天，实在想不到什么了，在看似山穷水尽的时候，和同事 Allen的聊天提醒了我。
“在求交点的两重循环代码中，你是把三角形放在外循环，还是把射线放在外循环？“
“我把射线放在外循环，有什么区别吗？”
“有区别，区别很大，要把三角形放在外循环“
“为什么？“
“局部性原理“
我一拍脑袋，顿时豁然开朗。对啊，程序中最耗时一部分就是频繁的通过做向量叉积计算三角形的法向量，把三角形放在外循环，就可以在内循环中复用已经计算好的三角形法向量。原先是先计算一条射线和所有三角形的交点，现在改为先计算一个三角形和所有射线的交点。这个简单的改变使程序运行时间几乎降了一半。其实就是基于操作系统的“局部性原理”，简而言之，就是说当前或者说最近被访问的指令和数据很可能被再次访问。优化的关键就是找出这些需要大量计算并且频繁访问的内容并尽可能复用结果。这个改变似乎开阔了我的思维，又有了优化的思路，比如在点积和叉积这些相对运算量较大的计算之前加一些简单的预算，能尽早跳出循环尽早跳出，避免做任何无意义的运算，尽量推迟、拆分、复用运算量较大的计算等等。改了一天程序，代码似乎变得难看了，但执行效率确有显著提高，下班前的测试结果使我兴奋异常---12秒多。我也接近最快的速度了。
第四天： SSE,是加速器还是减速带？
一大早还来不及沉浸在12秒的喜悦中，就被同事Jack当头一棒。2秒---一个恐怖的新纪录诞生了，由Jack创造。我的第一感觉是：太假了吧，是不是算错了。现实是Jack确实只用了2秒多完成那个最大的输入样本，而且结果正确。秘诀是他用SSE3替换了原本的C++代码。顿时，同事们都开始了解研究 SSE指令集。SSE是什么？其实我也是当时才知道。简单的说，SSE是Intel的SIMD扩展指令集，专门用于多媒体特别是图形学中大规模的矩阵运算和浮点运算。该指令较之高级语言优势在于<1>直接操作4个128位的寄存器，大大减少了内存与寄存器之间的IO消耗<2>一个寄存器可存储4个浮点数，并在一条指令中同时做四次浮点运算。对其有基本了解之后，参加竞赛的同事们都开始用SSE指令替换原本的C++代码。但现实并不尽如人意。许多同事抱怨用了 SSE不但结果不对，速度还反而慢了许多。调查原因，原来使用该指令集对内存的字节对齐也有要求，分配的内存必须是16字节对齐，要用declspec(align(16))关键字声明分配的内存。注意到这些细枝末节之后我也写了一个 SSE3版本的程序，但是发现速度几乎没有什么提高，仅仅提高了2秒，达到了10秒。SSE,你究竟是加速器还是减速带？2秒神话是如何创造出来的？
第五天： 一剑穿四心，关键在并行
从网上得知，SSE如果使用得当可以使代码运行效率提速三到四倍。可是我未见效率的提升，瓶颈在哪里呢？厚道的Jack公布了他的秘密：<1>一剑穿四心，一次计算一个三角形和四条射线的交点，充分发挥SSE指令的并行运算能力。<2>其次要了解各个SSE指令本身的性能，能用逻辑操作完成的就不要用算术运算，尽量用CPU执行起来更简单方便的指令达到目的。<3>避免SSE指令和高级语言穿插使用，因为这样会造成内存与寄存器之间的IO消耗，抵消了SSE的优势<4>而且要将x,y,z分量分开存储在三个寄存器中，分别计算结果，实践证明这一点也非常重要。在Jack的启发下，我终于用SSE3写了一个自认为比较完美的版本，最后单线程运行时间是1.5秒，双线程0.9秒。
总结： 实践出真知
完全沉浸在这个竞赛中的5天，使我获益良多，不仅在技术层面上，更重要的是在思想上。
在当前以互联网为主导的信息化社会，如何快速、有效地获取知识并能够通过自己的思考进行信息筛选非常重要。就这次竞赛来说，同样用google, 有人能够很快找到最优或接近最优的方法，后面的优化工作就小些；有的则会走很多弯路，比如我， 开始就用了运算很慢的方法，后面才找到更好的方法。
关于学习，“博”与“专”孰优孰劣的话题一直争论不休。在这次竞赛中，我深刻体会到知识广博的重要性和各个领域的互通性，除了基本的软件理论，还需要了解几何中的三角形重心坐标表示法，需要懂得基本的矢量变换，需要了解SIMD硬件指令集和基本的硬件架构才能高效的利用SSE指令。其实针对这个竞赛题目，还有一种据说效率更高，对SSE支持更好的算法，用到了普鲁克坐标系。可惜我现在还没时间去学习和了解这个。 在17，18世纪的欧洲，数学家往往也是物理学家，甚至是天文学家。这也说明了领域的互通和知识广博的重要性。
最后一点很重要，对于一个不确定的想法，不要用自己的主观设想先入为主的轻易否定它，实践，才能知道这个想法是否正确。这个道理看似简单，竞赛中我却犯了多次先入为主的错误。比如我开始不相信用C的fread/fwrite替换C++的iftream/oftream会提高运行效率，最后证实，在大量的IO操作中，C函数的效率优势还是很明显的。还有，输出结果集我开始是用stl的vector存储，优化时我自己写了一个vector,每次分配固定长度的内存块（比如1024），然后通过指针将所有前一个内存块和后一个内存块相连。其实就是链表结构，每个链表节点中又包含一个数组，数组用完后增加一个新的链表节点。这样可以避免stl中在超过vector已分配长度时发生的内存再分配，拷贝和原内存块的销毁操作。原以为用这个自定义的vector运行效率会提高很多，最后证实没什么提高…总之，实践出真知，求知的过程就是学习、思考、实践、反馈、调整不断循环往复的过程。
