在没有临时数组的情况下评估 Fortran 的 SUM(..)？

Question

可读性方面，我觉得写

更可取

momentum = sum( [( calculateMomentum(elements(i)), i=1, size(elements,1) )] )

超过

momentum = 0.0d0
do i = 1, size(elements,1)
    momentum = momentum + calculateMomentum(elements(i))
end do

因为第一个版本的形式是定义momentum的值，而第二个版本对应的是更底层的累加指令。在具有更多参数和可能多个索引的实际代码中，差异变得更加明显。

但是，第一个版本分配了一个临时数组。作为一名人类程序员，我知道它可以被优化掉，所以我想知道 Fortran 是否提供了一种语法，它允许既不使用显式循环也不使用临时数组来计算总和。

更新

没有这样的构造。 看起来没有语法，正如我所要求的。参见。

这很重要，但没有我想象的那么重要。我使用矩阵乘法和多个变体制定了自己的基准 (pastebin, embedded)。

C(i,j) = C(i,j) + A(i,k) * B(k,j) 是最慢的，可能是由于每一步中不必要的数组访问。
C(i,j) = sum( [(A(i,k) * B(k,j), k = 1, N)] )` 比速度快 10-20% (1)、尽管临时数组。
tmp = tmp + A(i,k) * B(k,j)，即使用临时累加器变量，比 (2) 快 20%。
C = matmul(A,B) 是迄今为止最快的，与 (3) 相比，500x500 矩阵增加了 25 倍，2000x2000 增加到 50，而其他变体的相对速度保持不变大致相同。

底线：当任务无法用优化的库或内部函数轻松表达时，sum 变体具有可行的性能，并且只应优化掉，如果性能真的如此重要那部分代码。

Answer 1

是否分配或净分配临时数组取决于编译器中的优化。无论如何，堆栈分配几乎是免费的。复制值可能需要更长的时间。

如果编译器可以确保结果相同，它可能会优化不必要的步骤。但是，没有特殊的语法。 Fortran 通常会尽量远离实际实现，并在编译器上留下很多东西。

Answer 2

为了实验，我尝试了这段代码（计算 arr 的倒数之和）。

program main
    use iso_fortran_env, only: dp => real64
    implicit none
    real(dp) val
    real(dp), allocatable :: arr(:)
    integer num, loop, i, t1, t2, trate

    num = 10**8
    arr = [( i, i = 1, num )]   !! L1

    do loop = 1, 10
        call system_clock( t1 )

        val = sum( [( testfunc( arr(i) ), i = 1, num )] )   !! L2

        call system_clock( t2, trate )
        print *, "val = ", val, " in ", (t2 - t1) / real(trate), " (s)"
    enddo

contains

function testfunc( x ) result( ret )
    real(dp), intent(in) :: x
    real(dp) :: ret

    ret = 1.0_dp / x
end

end program

然后，在我的电脑上 (mac2012)，“gfortran-10 -O2 test.f90 && time ./a.out”给出

 val =    18.997896413852555       in    1.02999997      (s)
 val =    18.997896413852555       in    1.10099995      (s)
 val =    18.997896413852555       in    1.17600000      (s)
 ...

real    0m12.575s
user    0m8.142s
sys 0m4.387s

和“gfortran-10 -O3”给出

 val =    18.997896413852555       in   0.875000000      (s)
 val =    18.997896413852555       in   0.888000011      (s)
 val =    18.997896413852555       in   0.833000004      (s)
 ...

real    0m9.986s
user    0m5.738s
sys 0m4.210s

在这两种情况下，htop 命令显示已分配约 1.5 GB，如果行 L1 和 L2 使用临时数组（每行约 800 MB，分配约 0.3 秒），这可能是合理的。

因为没有用于创建“迭代器”的语法，所以我尝试使 testfunc() 成为基本的（或不纯的基本的）。这里唯一的区别是标有 <--.

的行

program main
    use iso_fortran_env, only: dp => real64
    implicit none
    real(dp) val
    real(dp), allocatable :: arr(:)
    integer num, loop, i, t1, t2, trate

    num = 10**8

    arr = [( i, i = 1, num )]

    do loop = 1, 10
        call system_clock( t1 )

        val = sum( testfunc( arr ) )   !<--

        call system_clock( t2, trate )
        print *, "val = ", val, " in ", (t2 - t1) / real(trate), " (s)"
    enddo

contains

impure elemental &   !<--
function testfunc( x ) result( ret )
    real(dp), intent(in) :: x
    real(dp) :: ret

    ret = 1.0_dp / x
end

end program

然后，“gfortran-10 -O2”给出

 val =    18.997896413852555       in   0.437000006      (s)
 val =    18.997896413852555       in   0.453999996      (s)
 val =    18.997896413852555       in   0.437999994      (s)
 ...
real    0m5.946s
user    0m5.069s
sys 0m0.842s

和“gfortran-10 -O3”给出

 val =    18.997896413852555       in   0.225999996      (s)
 val =    18.997896413852555       in   0.252000004      (s)
 val =    18.997896413852555       in   0.246999994      (s)
 ...
real    0m3.909s
user    0m3.009s
sys 0m0.867s

htop 命令显示 ~800 MB，因此似乎只分配了 arr。

为了比较，以下代码使用显式 do-loop 计算 val（使用 testfunc() 的标量版本）

val = 0
do i = 1, num
    val = val + testfunc( arr(i) )
enddo

与上面 elemental + testfunc(arr) 的第二个代码给出了相同的时间（在我的 mac 上）。

在上述所有代码中，-Ofast 选项导致运行时错误 (Illegal instruction)。但这是由于行 L1 (arr = [( i, i = 1, num )])。如果我预先分配 arr 并用显式循环填充它，-Ofast 也可以正常工作（在这种情况下与 -O3 的时间几乎相同）。

Answer 3

如果你制作 calculateMomentum 和基本函数，那么它既可以用于标量值，也可以用于数组

例如：

    elemental function calculateMomentum(obj, v) result(p)
    class(body), intent(in) :: obj
    real, intent(in) :: v
    real :: p
        p = obj%mass * v
    end function

您可以将以上内容应用于 v

的数组

integer::i
type(body) :: ball
real, allocatable :: v(:), p(:)
real :: tot_p
allocate(v(10))
v = [ (10+i, i=1, 10) ]
p = calculateMomentum(ball, v)
tot_p = sum(p)

有一个中间数组来保存值是有利的，因为它使数据保持在附近（可能在 cache-line 内）并且 sum() 函数将尽可能快。

设想类型body例如如下

    type body
        real :: mass
    end type

在没有临时数组的情况下评估 Fortran 的 SUM(..)？

Evaluation of Fortran's SUM(..) without temporary array?

fortran

更新