2

我正在将 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 代码的循环顺序，或创建临时 2D 矩阵和 1D 数组以在使用 SUM 时“减少”维度，但显式 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
4

3 回答 3

0

首先，欢迎来到 StackOverflow。请参加游览！我们期望一个最小、完整和可验证的示例是有原因的，因为我们查看您的代码并且只能猜测可能是什么情况，这对社区没有太大帮助。

我希望以下建议可以帮助您弄清楚发生了什么。

使用 size() 函数并打印 Fortran 认为的尺寸大小以及打印 nx、ny 和 nz。据我们所知，数组被声明为大于 nx、ny 和 nz，并且这些变量是根据数据集设置的。Fortran 不一定将数组初始化为零，具体取决于它是静态数组还是可分配数组。

您还可以尝试在 sum 函数中指定数组范围：

r = Sum(foo(1:nx,1:ny,1:nz))

如果这样做，至少我们知道 sum 函数正在处理循环循环的完全相同的 foo 切片。

如果是这种情况，即使代码没有任何“错误”，您也会得到错误的答案。这就是为什么给出最小的、完整的和可验证的例子特别重要。

于 2019-10-06T15:59:05.597 回答
0

The sum intrinsic function returns a processor-dependant approximation to the sum of the elements of the array argument. This is not the same thing as adding sequentially all elements.

It is simple to find an array x where

summation = x(1) + x(2) + x(3)

(performed strictly left to right) is not the best approximation for the sum treating the values as "mathematical reals" rather than floating point numbers.

As a concrete example to look at the nature of the approximation with ifort, we can look at the following program. We need to enable optimizations here to see effects; the importance of order of summation is apparent even with optimizations disabled (with -O0 or -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

If adding up in strict order we get 1., seeing that anything smaller in magnitude than epsilon(0.) doesn't affect the sum.

You can experiment with the size of the array and order of its elements, the scaling of the small numbers and the ifort floating point compilation options (such as -fp-model strict, -mieee-fp, -pc32). You can also try to find an example like the above using double precision instead of default real.

于 2019-10-06T16:31:42.873 回答
0

我现在可以看到差异。这些是将小数加到大和的典型舍入误差。处理器可以使用它想要的任何求和顺序。没有“正确”的顺序。您不能真的说原始循环做出了“正确”的答案，而其他循环则没有。

您可以做的是使用double precision。在极端情况下，有像 Kahan summation 这样的技巧，但很少需要它。

将小数加到大和中是不精确的，尤其是在单精度中。您的结果中仍有四位有效数字。

通常不DIM=使用在某些特殊情况下使用的参数。

如果要对 的所有元素求和foo，请仅使用

s0 = SUM(foo)

足够了。

什么

s0 = SUM(SUM(SUM(foo, DIM=3), DIM=2), DIM=1)

所做的是它会创建一个临时二维数组，其中每个元素是 z 维中相应行的总和，然后是一个一维数组，每个元素是二维数组最后一维的总和，最后是该一维的总和大批。如果做得好，最终的结果将是相同的，但它会吃掉很多 CPU 周期。

于 2019-10-06T16:18:55.330 回答