计算结果差异产生来源

芯片架构差异

从底层架构来看，鲲鹏处理器和x86的芯片架构不同，如x86设计有40位、80位的浮点计算单元，在鲲鹏处理器上并没有这些计算单元。如果x86平台有使用到这些计算单元，就会导致两个平台的计算结果不一致。目前可以在鲲鹏平台上通过软件模拟的方式来解决这类问题引起的差异。例如icc的数学库使用了powr8i4函数，该函数调用了x87的80位指令，即用80位来表示浮点数。在鲲鹏平台上可以使用支持多种精度计算的开源MPFR（Multiple Precision Floating-Point Reliable）库来实现80位浮点数的pow运算。

FTZ（Flush-to-Zero）和DAZ（Denormals-are-Zero）：
FTZ和DAZ是专门用来处理非规格浮点数的标准，非规格浮点数定义解释见浮点数简介。FTZ表示浮点运算结果输出为非规格浮点数时转化为0，DAZ表示输入为非规格浮点数时转化为0。FTZ转化动作由硬件完成，能够加快非规范浮点数的计算速度。有如下示例代码：
```
float a = 1.0e-5 * 1.0e-38;
cout << "a = " << a << endl;
```
在不使能FTZ功能下，输出为：
```
a = 9.94922e-44
```
使能FTZ功能下，gcc增加-ffast-math编译选项即使能FTZ，输出为：
```
a = 0
```
查看增加-ffast-math编译选项后的汇编代码（鲲鹏平台上的代码），多了如下汇编指令：
```
0000000000400770 <set_fast_math>:
  400770:       52a02000        mov     w0, #0x1000000                  // #16777216
  400774:       d51b4400        msr     fpcr, x0
  400778:       d65f03c0        ret
  40077c:       00000000        .inst   0x00000000 ; undefined
```
通过赋值FPCR（Floating-Point Control Register）寄存器来控制FTZ的使能，鲲鹏和x86都支持FTZ功能，但由于两者的硬件架构不一样，存在一些细微的差别，鲲鹏和x86对于非规格浮点数执行FTZ和舍入的先后顺序不一样。有下面代码：
```
float a = 1.09628081709e-33;
float b = 1.07225660031e-05;
float c = a * b;
cout  << setprecision(12) << "c = " <<  c << endl;
```
鲲鹏机器下使能FTZ，g++ -ffast-math，输出结果为：
```
c = 0
```
x86机器下使能FTZ，g++ -ffast-math，输出结果为：
```
c = 1.17549435082e-38
```
1.17549435082e-38（0x800000）是IEEE 754标准中单精度能够表示的最接近0的浮点数，科学计算中1.09628081709e-33 x 1.07225660031e-05 = 1.17549434192e-38，计算结果比c值小，在x86机器上舍入到最接近的浮点数1.17549435082e-38（0x800000）。鲲鹏和x86执行动作对比如下。
- 鲲鹏执行顺序
  1. 执行乘法计算，a * b = 1.17549434192E-38。
  2. 执行FTZ，1.17549434192e-38为非规格浮点数，转化为0。
  3. 执行舍入，0舍入到最接近的浮点数，结果还是0。
  4. 执行赋值，将0赋值给变量c。
- x86执行顺序
  1. 执行乘法计算，a * b = 1.17549434192E-38。
  2. 执行舍入，1.17549434192e-38舍入到最接近的浮点数1.17549435082e-38（0x800000）。
  3. 执行FTZ，1.17549435082e-38（0x800000）不是非规格浮点数，不会转化为0。
  4. 执行赋值，将1.17549435082e-38（0x800000）赋值给变量c。
浮点数转换溢出：
浮点数转换为整数类型时，通常会保留整数部分，截取掉小数部分。当浮点数太大，整数类型无法表示浮点数的整数部分时，鲲鹏和x86的行为是不一致的，鲲鹏对于上溢或下溢的行为处理非常清晰和简单，即保留整数部分能表示的最大值或最小值。x86则定义了一个不确定数值（indefinite integer value）来表示溢出值。双精度转换的不确定数值为：0x8000 0000 0000 0000，单精度转换的不确定数值为：0x8000 0000。有下面示例代码：
```
double a = 1.0e20;
long b = static_cast<long>(a);
std::cout << "b = "  << b << std::endl;
```
鲲鹏平台下输出结果为：
```
b = 9223372036854775807
```
x86平台下输出结果为：
```
b = -9223372036854775807
```
上面示例中1.0e20超出long能够表示的最大值，鲲鹏平台上把long能够表示的最大整数赋值给b = 9223372036854775807(0x7FFF FFFF FFFF FFFF)。x86平台上把双精度转换的不确定数值赋值给b = -9223372036854775807(0x8000 0000 0000 0000)。

部分浮点数转换溢出时鲲鹏和x86的结果对比请参考《TaiShan服务器代码移植参考》案例“双精度浮点型转整型时数据溢出，与x86平台表现不一致”的相关内容。

编译器差异

编译选项-fp-model：
icc编译器（Intel C++ Compiler）的编译选项-fp-model是用来控制浮点数是否值安全（value-safe），即编译器是否进行在数学逻辑上相等的运算转换，当浮点运算不是值安全时可能会产生计算结果差异。-fp-model默认的选项是fast=1，在-O2及以上优化等级上默认开启。可以指定-fp-model=precise来保证浮点运算是值安全的。
浮点数无法准确表示实数，鲲鹏和x86的默认舍入方式都采用舍入到最接近的实数。舍入带来精度损失，会导致计算机的运算不一定符合数学逻辑上的运算。例如A + B + C 不一定等于 A + ( B + C )，如下面代码所示：
```
float a = -0.1;
float b = 0.1;
float c = -0.00000000149012;
float sum1 = a + b + c;
float sum2 = a + (b + c);
std::cout << "sum1 = " << sum1 << std::endl;
std::cout << "sum2 = " << sum2 << std::endl;
```
输出结果为：
```
sum1 = -1.49012e-09
sum2 = 0
```
按照数学运算结果sum1应该等于sum2，问题出现在(b + c)计算上，c(-0.00000000149012)约等于b(0.1)的ULP的五分之一，ULP定义解释见浮点数简介，按照舍入到最近的数的方式，c将会被舍去，(b + c) = b，因此sum2 = a + b = -0.1 + 0.1= 0。

编译器会对代码进行性能优化，有时会隐式打乱计算顺序。而不同的编译器有不同优化策略，会导致在一样的代码下，使用不同的编译器运行结果也会不同。有下面示例代码：
```
#include <iostream>
float a[4 * 4] = {0};
float sum = 0;
void main ()
{
    a[0] = 0.1;
    a[1] = -0.1;
    a[2] = 0.00000000149012;
    for (int i = 0; i < 4 * 4; i++) {
        sum += a[i];
    }
    std::cout << "sum = " << sum << std::endl;
}
```
在icc -O3条件下，输出结果为：
```
sum = 0
```
在g++ -O3条件下，输出结果为：
```
sum = 1.49012e-09
```
icc的计算结果跟上个示例中a + (b + c)结果一样。因为在-O3优化等级下，icc默认开启-fp-model=fast，对a[i]进行了向量化优化时，进行了符合数学逻辑运算的转换，导致计算没有按照a[0] + a[1] + a[2]...顺序进行，而gcc编译器则没有进行向量化优化，保持a[0] + a[1] + a[2]...的计算顺序，所以两者结果出现不一致的现象。

在icc -O3 -fp-model=precise条件下，输出结果为：
```
sum = 1.49012e-09
```
icc编译器在-fp-model=precise使能情况下，会禁止包括但不限于下面的数学逻辑相等的运算转换：

a + b + c → a + (b + c) 、x * 0 → 0 、x + 0 → 0 、A / B → A * (1 / B)
立即数：
立即数在代码中是显示表示的数，对于未定义类型的立即数，其加载到内存的值在IEEE 754中未作具体的定义，因此不同的编译器对立即数加载有不同实现时，有时会导致立即数在内存的数值不一致，有下面示例代码：
```
program test
real (8) :: x
x = 0.002362E12
write(*, '(Z16)') x
end
```
使用gcc编译，输出为：
```
41E1992840000000
```
使用icc编译，输出为：
```
41E1992860000000
```
指定立即数类型为双精度浮点数就不会出现差异：
```
program test
real (8) :: x
x = 0.002362E12_8 ! 指定类型为双精度
write(*, '(Z16)') x
end
```
使用gcc编译，输出为：
```
41E1992850000000
```
使用icc编译，输出为：
```
41E1992850000000
```

数学库差异

如果函数在任何情况下产生的误差总是小于0.5ulp，那么说明这个函数总是能够返回最接近精确结果的近似值。但是目前没有统一的标准去约束数学库中log、sin、pow等函数的输出，因此不同的数学库中的函数实现可能具有不一样的舍入模式和精度计算。有下面示例代码：

program test
real::x, y, z
x = 0.699129045
y = 1.4
z = x ** y
print*, z
end

使用系统数学库libm.so.6编译，输出结果为：

0.605871201

使用鲲鹏数学库libkm.so编译，输出结果为：

0.605871141

应用代码引起的差异

WRF网格划分引起的差异，有如下代码：
```
ht_loc(its-1, j) = max(ht_loc(its-1, j), ht_shad(its-1, j))
ht_loc(ite-1, j) = max(ht_loc(ite-1, j), ht_shad(ite-1, j))
ht_loc(i, jts-1) = max(ht_loc(i, jts-1), ht_shad(i, jts-1))
ht_loc(i, jte-1) = max(ht_loc(i, jte-1), ht_shad(i, jte-1))
```
代码中的its、ite、jts、jte变量表示网格区域范围，其划分范围跟使用的进程数有关，进程数越多时，每个进程分配的网格区域范围越小，反之分配的网格区域范围越大。当更改测试规模时（测试节点数、测试线程数），its、ite、jts、jte变量值发生改变，变量ht_loc也一同改变，经过多次迭代运算后，整个程序的输出结果也发生改变。
代码中的未定义行为，如数组越界bug，有如下代码：
```
real :: SLA_TABLE(27)
// i为外部输入
masslai =  max(SLA_TABLE(i), 1.0)
```
代码中SLA_TABLE的数组大小为27，变量i是外部输入，当i大于27时，数组发生越界，SLA_TABLE(i)的值不可确定，masslai变量的值也不可确定，导致最终输出结果也不确定。