gromacs, quake III和vrsqrtps

时间:2023-03-09 20:20:25
gromacs, quake III和vrsqrtps

看标题大家可能觉得三个词汇风马牛不相及,第一个是解蛋白质分子动力学的软件,第二个是上三代宅男最爱雷神之锤,第三个则是一个存在于IntelSSE及AVX中的一个指令,他的作用是快速求平方根的倒数。

起因是这样子的。某天闲着没事,跑去benchmarksgame.alioth.debian.org上看到了万年被压在fortran身体下蹂躏的c++居然翻身了。最不可思议的是,在fortran长项上的多体运算nbody居然被c++拉了一大截性能下来(2倍,5千万步,fortran用时19秒,C++用时9.59秒)。本着fortran居然能被c++搞倒根本是不可能的事情的信念,我开始出发了。

好在他们提供了nbody的f90代码。装好ifort就开始优化,初一眼看过去就好多问题,“怪不得那么慢”,看来现在fortran培训太少了以至于部分不合格码农直接拿java代码改改就过来用了。内存不对齐,行列优先搞混,没有写成向量化,这类错误比比皆是。改了改,差不多跑快了10秒,感觉差不多了,把c++的代码拉下来一编译,瞬间崩溃,在本机上只跑了6秒!

接下去是苦逼的两天,我绞尽脑汁改代码,向量化,发现部分代码无法向量化,然后手工指定,变更慢了……终于,在现实面前,终于低下头,想想是不是c++里用了什么新技术导致比fortran快呢?或许人家也有新的算法呢?(可喜可贺)

跑去看c++的代码部分。就发现他将多体问题中常见的求出距离求倒数r^(-1)=(1/sqrt(x*x+y*y+z*z))的那一步用了很简明的方法写了一下

distance=_mm_cvtps_pd(_mm_rsqrt_ps(_mm_cvtpd_ps(dsquared)));

这里给不是很了解这段代码的人解释一下,dsquared是一个双精度浮点数,_mm_cvtpd_ps是把一个双精度浮点数变成了单精度浮点数。_mm_cvtps_pd是把一个单精度浮点数变成双精度浮点数。_mm_rsqrt_ps就难倒我了,不懂c++的人表示很不解:这是什么东西啊?拿去搜了一把发现这其实是调用SSE指令集vrsqrtps的一个函数,好高端!居然能在C++里直接实现SSE/AVX指令集的调用。再一搜索,差点笑喷出来,搞了半天敢情这个算法就是QuakeIII中开用的快速求平方根倒数的那个算法。也就是传说中的超快算法,利用泰勒展开来求解。所以在intel官网上,关于这一指令的介绍是这么写的“(V)RSQRT[P/S]SComputeapproximatereciprocalofsquarerootofpacked/scalarsingleprecision”

关于这个fastinversesquareroot的具体步骤可以看wiki:http://en.wikipedia.org/wiki/Fast_inverse_square_root

但很遗憾的是,这只是一种近似算法。在大多数情况下一次迭代能收敛到单精度准确值。但二次迭代就变慢了,所以很多人不使用二次迭代。雷神之锤里用的也是一次迭代,而对于一个3d游戏来说这个精度足够了。在浏览之余,我看到了某个小论坛上表明GROMACS其实也是使用这个算法的,于是就去搜了一下,在最新版本4.6.3中也发现了类似代码(以下只是专门为avx指令集编写的头文件中的一个,有兴趣的人可以自己在gromacs源代码目录中用grep-IRn_mm_rsqrt_ps*看一下)

include/gmx_math_x86_avx_128_fma_double.h:

/* 1.0/sqrt(x) */
static gmx_inline __m128d
gmx_mm_invsqrt_pd(__m128d x)
{
const __m128d half = _mm_set1_pd(0.5);
const __m128d three = _mm_set1_pd(3.0); /* Lookup instruction only exists in single precision, convert back and forth... */
__m128d lu = _mm_cvtps_pd(_mm_rsqrt_ps( _mm_cvtpd_ps(x))); lu = _mm_mul_pd(_mm_mul_pd(half, lu), _mm_nmacc_pd(_mm_mul_pd(lu, lu), x, three));
return _mm_mul_pd(_mm_mul_pd(half, lu), _mm_nmacc_pd(_mm_mul_pd(lu, lu), x, three));
}

不同之处在于,这个算法是完整的fastinversesquareroot算法(做了一步迭代以保证收敛到单精度),而nbody中的那个fastinverseSquareroot算法则没有做这一步迭代。这里要提醒诸位注意。这一函数对于双精度程序也会降低精度到单精度。由于_mm_cvtpd_ps已经把一个双精度数值降低到单精度了,所以_mm_cvtps_pd拿到的是一个单精度值,而对于单精度数值变成双精度数值,计算机的做法就是剩余位数用随机数填充。所以无论是gromacs编译时如何指定单双精度。拿到距离矩阵已经是单精度了。

接下去做一个比较

双精度1/sqrt

~$./nbody2f50000000
初始能量-0.1690751638285244717874178377
终态能量-0.1690599067877010530658310472

单精度1/sqrt
~$./nbodyc50000000
初始能量-0.1690751638285244717874178377
终态能量-0.1690599067881611849983869433

加黑部分就是误差了。而且这个体系非常小,只有5体运动,而对于一般性蛋白质模拟,可能牵涉到几百万个原子,这个误差积累就相当大了

然后,这里测试了5亿个数,从0开始步长0.0001的1/sqrt测试(也就是0.0001,0.0002,0.0003.....到50000)。我发现在不迭代时,最大误差能达过8.5e-3。当计算3e-4的1/sqrt时,误差能达到8.5e-3,平均误差则在3.57e-7,IEEE允许的单精度误差是1.17e-7,而双精度允许误差则是2.22e-16。所以说不包含一次迭代的误差已经超过单精度允许误差了。

而包含一次迭代的算法,最大误差在2.7e-13,高于双精度允许误差但低于单精度允许误差,有问题的是2.1e-3的1/sqrt(2.1e-3)。平均误差8.67e-19。

从误差分布来,1/sqrt(x)的快速算法中,x与误差是成反比的,x越靠近0,误差越大,x越大误差则越小。

到这里。谜底揭晓了,alioth上那个速度超快的c++程序和gromacs都是通过降低精度来达到高速的,fortran的王者地位依然无法动摇(偷笑)。

Reference

http://software.intel.com/en-us/articles/introduction-to-intel-advanced-vector-extensions

ftp://ftp.gromacs.org/pub/manual/manual-4.5.6.pdf AppendixB.3,我猜没人会看到附录B.3的……

http://scicomp.stackexchange.com/questions/2168/what-is-the-computational-cost-of-sqrtx-in-standard-libraries