我正在将 f77 代码转换为 f90 代码，部分代码需要对 3d 矩阵的元素求和。在 f77 中，这是通过使用 3 个循环(外部、中间、内部索引(来实现的。我决定使用 f90 内禀和(3 倍(来实现这一点，令我惊讶的是答案不同。我正在使用ifort编译器，已打开调试，检查边界，未优化
这是 f77 样式的代码
r1 = 0.0
do k=1,nz
do j=1,ny
do i=1,nx
r1 = r1 + foo(i,j,k)
end do
end do
end do
这是 f90 代码
r = SUM(SUM(SUM(foo, DIM=3), DIM=2), DIM=1)
我尝试了各种变化，例如交换 f77 代码的循环顺序，或者在使用 SUM 时创建临时 2D 矩阵和 1D 数组以"减少"维度，但显式 f77 样式循环总是给出与 f90+ SUM 函数不同的答案。
我将不胜感激任何有助于理解差异的建议。
顺便说一下，这是使用一个串行处理器。
编辑 12：13 下午以显示完整示例
! ifort -check bounds -extend-source 132 -g -traceback -debug inline-debug-info -mkl -o verify verify.f90
! ./verify
program verify
implicit none
integer :: nx,ny,nz
parameter(nx=131,ny=131,nz=131)
integer :: i,j,k
real :: foo(nx,ny,nz)
real :: r0,r1,r2
real :: s0,s1,s2
real :: r2Dfooxy(nx,ny),r1Dfoox(nx)
call random_seed
call random_number(foo)
r0 = 0.0
do k=1,nz
do j=1,ny
do i=1,nx
r0 = r0 + foo(i,j,k)
end do
end do
end do
r1 = 0.0
do i=1,nx
do j=1,ny
do k=1,nz
r1 = r1 + foo(i,j,k)
end do
end do
end do
r2 = 0.0
do j=1,ny
do i=1,nx
do k=1,nz
r2 = r2 + foo(i,j,k)
end do
end do
end do
!*************************
s0 = 0.0
s0 = SUM(SUM(SUM(foo, DIM=3), DIM=2), DIM=1)
s1 = 0.0
r2Dfooxy = SUM(foo, DIM = 3)
r1Dfoox = SUM(r2Dfooxy, DIM = 2)
s1 = SUM(r1Dfoox)
s2 = SUM(foo)
!*************************
print *,'nx,ny,nz = ',nx,ny,nz
print *,'size(foo) = ',size(foo)
write(*,'(A,4(ES15.8))') 'r0,r1,r2 = ',r0,r1,r2
write(*,'(A,3(ES15.8))') 'r0-r1,r0-r2,r1-r2 = ',r0-r1,r0-r2,r1-r2
write(*,'(A,4(ES15.8))') 's0,s1,s2 = ',s0,s1,s2
write(*,'(A,3(ES15.8))') 's0-s1,s0-s2,s1-s2 = ',s0-s1,s0-s2,s1-s2
write(*,'(A,3(ES15.8))') 'r0-s1,r1-s1,r2-s1 = ',r0-s1,r1-s1,r2-s1
stop
end
!**********************************************
sample output
nx,ny,nz = 131 131 131
size(foo) = 2248091
r0,r1,r2 = 1.12398225E+06 1.12399525E+06 1.12397238E+06
r0-r1,r0-r2,r1-r2 = -1.30000000E+01 9.87500000E+00 2.28750000E+01
s0,s1,s2 = 1.12397975E+06 1.12397975E+06 1.12398225E+06
s0-s1,s0-s2,s1-s2 = 0.00000000E+00-2.50000000E+00-2.50000000E+00
r0-s1,r1-s1,r2-s1 = 2.50000000E+00 1.55000000E+01-7.37500000E+00
首先，欢迎来到 StackOverflow。请参加参观！我们期待一个最小、完整和可验证的示例是有原因的，因为我们查看您的代码，只能猜测可能的情况，这对社区没有太大帮助。
我希望以下建议可以帮助您弄清楚发生了什么。
使用 size(( 函数并打印 Fortran 认为的尺寸以及打印 nx、ny 和 nz。据我们所知，数组被声明为大于 nx、ny 和 nz，这些变量是根据数据集设置的。Fortran 不一定将数组初始化为零，具体取决于它是静态数组还是可分配数组。
您也可以尝试在 sum 函数中指定数组扩展数据块：
r = Sum(foo(1:nx,1:ny,1:nz))
如果这样做，至少我们知道 sum 函数正在处理循环遍历的完全相同的 foo 切片。
如果是这种情况，即使代码没有任何"错误"，您也会得到错误的答案。这就是为什么给出最小、完整和可验证的例子特别重要的原因。
我现在可以看到差异。这些是将小数字相加到大和的典型舍入误差。允许处理器使用它想要的任何求和顺序。没有"正确"的顺序。你不能真的说原始循环做出"正确"的答案，而其他循环则没有。
您可以做的是使用双精度。在极端情况下，有像卡汉求和这样的技巧，但很少需要它。
将小数加到大和是不精确的，在单精度中尤其如此。结果中仍有四位有效数字。
通常不使用在某些特殊情况下使用的DIM=参数。
如果要对foo的所有元素求和，只需使用
s0 = SUM(foo)
这就够了。
什么
s0 = SUM(SUM(SUM(foo, DIM=3), DIM=2), DIM=1)
这样做是，它将创建一个临时 2D 数组，其中每个元素是 z 维度中相应行的总和，然后是一个 1D 数组，每个元素是 2D 数组最后一个维度的总和，最后是该 1D 数组的总和。如果做得好，最终结果会是一样的，但它会消耗大量的CPU周期。
sum内部函数返回与处理器相关的近似值，以接近数组参数的元素之和。 这与按顺序添加所有元素不同。
找到一个数组很简单x其中
summation = x(1) + x(2) + x(3)
(严格从左到右执行(不是将值视为"数学实数"而不是浮点数的总和的最佳近似值。
作为查看 ifort 近似性质的具体示例，我们可以查看以下程序。 我们需要在此处启用优化才能看到效果;即使禁用了优化(使用-O0或-debug(，求和顺序的重要性也很明显。
implicit none
integer i
real x(50)
real total
x = [1.,(EPSILON(0.)/2, i=1, SIZE(x)-1)]
total = 0
do i=1, SIZE(x)
total = total+x(i)
print '(4F17.14)', total, SUM(x(:i)), SUM(DBLE(x(:i))), REAL(SUM(DBLE(x(:i))))
end do
end program
如果按照严格的顺序加起来，我们得到1.，看到任何小于epsilon(0.)的量级都不会影响总和。
您可以尝试数组的大小及其元素的顺序、小数的缩放和 ifort 浮点编译选项(如-fp-model strict、-mieee-fp、-pc32(。 您也可以尝试使用双精度而不是默认的 real 来查找类似上面的示例。