主页 › 博客 ›

关于科学计算中的数值误差问题

作者：甘驰 (Intel) (16 篇文章) 日期：一月 26, 2010 在 3:51 下午

最近看到有开发HPC应用的程序员反映换了编译器后，新程序的结果不正确。我不禁想起20年前在16位IBM PC AT/XT台式机上开发CAD/CAM程序, 程序在计算边界状况时，如两圆柱的轴夹角为1分求交线，结果经常匪夷所思，都是精度惹得祸。进入64bit时代同样的问题仍然存在，单次运算的精度提高了，但HPC程序中的大量的迭代计算，累积了误差，导致结果不同。我读了Intel同事的文章（http://software.intel.com/en-us/articles/consistency-of-floating-point-results-using-the-intel-compiler），略有所得，与大家分享。

人们要求结果的复现性有以下几方面，
-同一程序，同一数据集，执行若干次的结果相同（近）
-同一代码，不同编译选项产生的程序，相同的数据集，执行的结果相同（近）
-同一代码，不同编译器产生的程序，相同的数据集，执行的结果相同（近）
-同一代码，相同的数据集，在不同架构的CPU上执行的结果相同（近）
最好运行速度也要快。这些要求有时是相互矛盾的。让我们先分析计算结果有差别的原因，注意这里我用了差别，而不是误差。

首先，通常十进制的浮点数用二进制数表示时就有误差，便有了‘原罪’。

其二，计算机运算是不符合结合律的，如A+B+C不一定等于A+(B+C)。当A= -B且C为极小值时，A+B+C=C，但A+(B+C)=0，极小值C和B相加时由于精度限制，被忽略了或者说表示B的有效位不能包含C的值；但A+B=0,再加C后，结果为C。当你多线程或多进程化程序运行时，原来的数据被切成几分，分给不同的进程/线程执行，于是原先的次序被打破，用MPI或OpenMP的程序都有引入的这类潜在的结果差别。

其三，X87 FPU是以80位（1位符号，15位指数，64位数值）、扩展双精度浮点数运算的。SIMD执行单元可能不是（查遍手册，未见答案，至少其数据是以64位存放的）。

其四，四舍五入，(a) 1.0001 0000 1000 0011 1001 0111E2 101以单精度数表示时需要四舍五入，因为其数值有25位，超过单精度浮点数的24位。根据系统设置，有四种模式可选，
-Round to nearest (even) --- 选精度误差最小的
-Round down (toward −∞) -- 选精度误差最小的, 但不大于原数
-Round up (toward +∞) --选精度误差最小的, 但不小于原数
-Round toward zero (Truncate) --选精度误差最小的, 但其绝对值不大于原数

其五，Flush-to-zero 和 denormal-to-zero. 如设定flush-to-zero, 当浮点运算发生underflow时，返回0值而不是像IEEE754要求的，返回非正规数据（denormalized data），这主要是为了提升性能。Denormal-to-zero是将非正规数据操作数置零，也是为性能。

其六，不良的编程习惯。如变量不置初值等

我的同事将随后发如何巧用Intel编译器的选项避免结果差别，和如何发现导致结果差别的代码。

分类：并行计算, 英特尔® 软件网络 2.0