为什么Fortran内部函数“扩散”通常比显式迭代慢

Question

我使用地球物理模型，常见的情况是需要将3D数据与3D数据相乘，相加等。以下是一个例子。

module benchmarks
  implicit none
  integer, parameter :: n=500
  integer :: k
  real :: d2(n,n)
  real :: d3(n,n,n)
  contains
  ! Iteration
  subroutine benchmark_a(res)
    real, intent(out) :: res(n,n,n)
    do k = 1, size(d3,3)
      res(:,:,k) = d2*d3(:,:,k)
    end do
  end subroutine
  ! Spread
  subroutine benchmark_b(res)
    real, intent(out) :: res(n,n,n)
    res = d3*spread(d2, 3, size(d3,3))
  end subroutine
end module

program main
  use benchmarks
  real :: t, tarray(2)
  real :: res(n,n,n)
  call random_number(d2)
  call random_number(d3)
  ! Iteration
  call dtime(tarray, t)
  call benchmark_a(res)
  call dtime(tarray, t)
  write(*,*) 'Iteration', t
  ! Spread
  call dtime(tarray, t)
  call benchmark_b(res)
  call dtime(tarray, t)
  write(*,*) 'Spread', t
end program

当我使用不同尺寸大小n运行时，我通常发现spread要慢得多;例如：

Spread   2.09942889
Iteration  0.458283991

有谁知道为什么spread接近而不是明确的for循环（我认为通常是不惜一切代价避免）是如此慢？

Answer 1

这里的基本答案是“它不是”。也许对于特定的编译器和特定情况，内在函数不像显式DO循环那样优化，但它不一定是那样。我使用ifort 19进行了测试，即使在默认优化级别，SPREAD内部和显式循环也生成类似的代码，当我更正程序以使用结果时，内在更快。

Iteration  0.2187500      0.1376885
 Spread  9.3750000E-02  0.1376885

我也会提醒（正如我在你的问题的评论中所做的那样）简单的基准程序通常不会衡量作者认为他们做了什么。您的原始和修订示例都表现出的最常见错误是，从未使用过测试工作的结果，因此足够聪明的编译器可以简单地消除整个操作。实际上，当我使用ifort 19构建两个测试用例时，编译器会完全删除所有工作，只留下时序代码。不用说，运行速度非常快。

  implicit none
  integer, parameter :: n=500
  integer :: k
  real :: d2(n,n)
  real :: d3(n,n,n)
  contains
  ! Iteration
  subroutine benchmark_a(res)
    real, intent(out) :: res(n,n,n)
    do k = 1, size(d3,3)
      res(:,:,k) = d2*d3(:,:,k)
    end do
  end subroutine
  ! Spread
  subroutine benchmark_b(res)
    real, intent(out) :: res(n,n,n)
    res = d3*spread(d2, 3, size(d3,3))
  end subroutine
end module

program main
  use benchmarks
  real :: tstart,tend
  real :: res(n,n,n)
  call random_number(d2)
  call random_number(d3)
  ! Iteration
  call cpu_time(tstart)
  call benchmark_a(res)
  call cpu_time(tend)
  write(*,*) 'Iteration', tend-tstart, res(10,10,10)
  ! Spread
  call cpu_time(tstart)
  call benchmark_b(res)
  call cpu_time(tend)
  write(*,*) 'Spread', tend-tstart, res(10,10,10)
end program```