Fortran 90数组双精度设置及硬编码赋值精度异常问题求助

问题背景 #

为了和C代码做数值对比，需要使用双精度数组。定义了如下双精度数组：

INTEGER :: n = 3
DOUBLE PRECISION, DIMENSION(n, 1) :: grid

随后通过硬编码赋值：

grid(1,1) = -2.85697001387280    
grid(2,1) = -1.35562617997427   
grid(3,1) =  0.00000000000000       

但调试时发现数组实际值为：

grid = [-2.8569700717926025, -1.3556262254714966, 0]

精度仅达到单精度水平。尝试用REAL*8、REAL(8)替换DOUBLE PRECISION后问题依旧。但通过计算生成的数组却能保持双精度，且和C代码数值一致：

INTEGER :: n = 3
REAL(8) :: mina = 0.5, maxa = 5.0 
REAL(8) :: lmina, step
REAL(8), DIMENSION(n,1) :: a

lmina = log(mina)
lmaxa = log(maxa)
step = (lmaxa - lmina) / (n - 1)
DO i=1,n
   a(i,1)=lmina+(i-1.0)*step
END DO

编译命令为：

gfortran -g -o myfile myfile.f90

原因分析 #

核心问题出在硬编码的字面量默认类型：Fortran中，不带后缀的小数字面量（比如-2.85697001387280）默认是单精度实数。当把单精度字面量赋值给双精度变量时，会先将单精度值转换为双精度，但单精度本身只有约7位有效数字，无法容纳输入的14位有效数字，转换过程中就会丢失精度，最终得到的双精度值其实是单精度值的“扩展”，自然达不到双精度的精度要求。

而计算生成数组的场景中，所有参与运算的变量都是双精度类型（REAL(8)），运算过程全程保持双精度，因此结果能达到双精度精度。

解决方法 #

给硬编码的小数字面量加上双精度后缀，让字面量本身成为双精度类型，这样赋值时就不会有精度损失。常用的双精度后缀是d0（大小写均可，D0也可以）：

grid(1,1) = -2.85697001387280d0    
grid(2,1) = -1.35562617997427d0   
grid(3,1) =  0.00000000000000d0       

这样赋值后，数组元素就能精确保留输入的双精度数值，和C代码的数值对比也能匹配。

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