同次运行循环内浮点运算一致性疑问——Intel Fortran相关

是的，Intel Fortran经典编译器的默认优化确实可能导致同一次运行、同一循环内，对相同输入的浮点运算产生不一致的结果，具体原因如下：

• 编译器的默认优化（如自动向量化、指令调度、寄存器重用等）会改变浮点运算的执行顺序或选用不同的指令集（如SSE/AVX）。不同指令集在中间结果存储精度上存在差异，比如部分运算可能临时使用80位x87寄存器存储中间值，另一部分则直接用64位双精度存储，这会导致相同输入的计算结果出现细微舍入差异。
• 循环迭代可能被编译器做并行化或向量化拆分处理，不同块的迭代运算路径可能被优化调整，进而出现舍入误差的不一致性。

针对你的情况，可以尝试以下解决办法：

• 开启Intel Fortran的-fp-model strict选项，强制编译器遵循严格浮点运算规则，减少优化带来的精度差异，但会牺牲部分性能。
• 对原式进行代数化简，减少舍入误差概率：
  原式：

  y = (p*q*x+ q*q/2)/(sqrt(p*p+0.25_dp*q*q+p*q*x)*q)
  

  化简后：

  y = (2*p*x + q) / (2*sqrt(p*p + p*q*x + 0.25_dp*q*q))
  

  该式减少了乘法次数，也避免了分母与分子中q的抵消误差。
• 在计算后添加clamp操作，规避微小误差带来的后续问题：

  y = max(min(y, 1.0_dp), -1.0_dp)
  

内容的提问来源于stack exchange，提问作者Jasper