numpy 怎么能比我的 Fortran 例程快得多？

Posted 2023-02-23

【中文标题】numpy 怎么能比我的 Fortran 例程快得多？

【英文标题】：How can numpy be so much faster than my Fortran routine?

【发布时间】：2016-02-16 20:47:20

【问题描述】：

我得到一个 512^3 数组，表示来自模拟的温度分布（用 Fortran 编写）。该数组存储在一个大小约为 1/2G 的二进制文件中。我需要知道这个数组的最小值、最大值和平均值，因为无论如何我很快就会需要理解 Fortran 代码，所以我决定试一试，并想出了以下非常简单的例程。

  integer gridsize,unit,j
  real mini,maxi
  double precision mean

  gridsize=512
  unit=40
  open(unit=unit,file='T.out',status='old',access='stream',&
       form='unformatted',action='read')
  read(unit=unit) tmp
  mini=tmp
  maxi=tmp
  mean=tmp
  do j=2,gridsize**3
      read(unit=unit) tmp
      if(tmp>maxi)then
          maxi=tmp
      elseif(tmp<mini)then
          mini=tmp
      end if
      mean=mean+tmp
  end do
  mean=mean/gridsize**3
  close(unit=unit)

在我使用的机器上，每个文件大约需要 25 秒。这让我觉得相当长，所以我继续在 Python 中执行以下操作：

    import numpy

    mmap=numpy.memmap('T.out',dtype='float32',mode='r',offset=4,\
                                  shape=(512,512,512),order='F')
    mini=numpy.amin(mmap)
    maxi=numpy.amax(mmap)
    mean=numpy.mean(mmap)

现在，我当然希望这会更快，但我真的被震撼了。在相同条件下只需不到一秒钟。平均值偏离了我的 Fortran 例程发现的平均值（我也使用 128 位浮点数运行，所以我更相信它）但仅在第 7 个有效数字左右。

numpy 怎么这么快？我的意思是您必须查看数组的每个条目才能找到这些值，对吗？我是否在我的 Fortran 例程中做了一些非常愚蠢的事情，以至于需要更长的时间？

编辑：

回答cmets中的问题：

是的，我还使用 32 位和 64 位浮点数运行了 Fortran 例程，但它对性能没有影响。 我使用了iso_fortran_env，它提供了 128 位浮点数。 使用 32 位浮点数虽然我的意思有点偏离，所以精度确实是个问题。 我以不同的顺序在不同的文件上运行了这两个例程，所以我猜在比较中缓存应该是公平的？ 我实际上尝试过打开MP，但要同时从不同位置的文件中读取。阅读了您的 cmets 并回答了这听起来现在真的很愚蠢，而且这也使例行程序花费了更长的时间。我可能会尝试一下数组操作，但也许这甚至没有必要。 文件大小实际上是 1/2G，这是一个错字，谢谢。 我现在将尝试数组实现。

编辑 2：

我实现了@Alexander Vogt 和@casey 在他们的答案中建议的内容，它和numpy 一样快，但现在我遇到了@Luaan 指出的精度问题。使用 32 位浮点数组，sum 计算的平均值为 20%。正在做

...
real,allocatable :: tmp (:,:,:)
double precision,allocatable :: tmp2(:,:,:)
...
tmp2=tmp
mean=sum(tmp2)/size(tmp)
...

解决了这个问题，但增加了计算时间（不是很多，但很明显）。 有没有更好的方法来解决这个问题？我找不到将文件中的单曲直接读取为双打的方法。 而numpy 又是如何避免这种情况的呢？

感谢到目前为止的所有帮助。

【问题讨论】：

您是否尝试过没有 128 位浮点数的 Fortran 例程？我不知道有任何硬件实际上支持这些，所以它们必须在软件中完成。

如果您尝试使用数组的 Fortran 版本（特别是使用一次读取而不是十亿次读取）会怎样？

您是否考虑过在 Fortran 中也使用数组运算符？然后，您可以尝试minval()、maxval() 和sum()？此外，您将 IO 与 Fortran 中的操作混合在一起，而不是在 Python 中 - 这不是一个公平的比较；-)

在对涉及大文件的内容进行基准测试时，请确保所有运行的缓存都相同。

另请注意，精度在 Fortran 中是一个相当大的问题，而且它是有代价的。即使你用你的 Fortran 代码解决了所有这些明显的问题，也很可能需要额外的精度，并且会导致显着的速度损失。

【参考方案1】：

您的 Fortran 实现有两个主要缺点：

您将 IO 和计算混合在一起（并逐个从文件条目中读取）。 您不使用向量/矩阵运算。

此实现确实执行与您相同的操作，并且在我的机器上快了 20 倍：

program test
  integer gridsize,unit
  real mini,maxi,mean
  real, allocatable :: tmp (:,:,:)

  gridsize=512
  unit=40

  allocate( tmp(gridsize, gridsize, gridsize))

  open(unit=unit,file='T.out',status='old',access='stream',&
       form='unformatted',action='read')
  read(unit=unit) tmp

  close(unit=unit)

  mini = minval(tmp)
  maxi = maxval(tmp)
  mean = sum(tmp)/gridsize**3
  print *, mini, maxi, mean

end program

想法是一口气将整个文件读入一个数组tmp。然后，我可以直接在数组上使用函数MAXVAL、MINVAL 和SUM。

---

对于准确性问题：只需使用双精度值并即时进行转换

mean = sum(real(tmp, kind=kind(1.d0)))/real(gridsize**3, kind=kind(1.d0))

只会略微增加计算时间。我尝试按元素和切片执行操作，但这只会增加默认优化级别所需的时间。

在-O3，逐元素加法的性能比数组运算好约 3 %。在我的机器上，双精度和单精度操作之间的差异小于 2% - 平均而言（单个运行偏差更大）。

---

这是一个使用 LAPACK 的非常快速的实现：

program test
  integer gridsize,unit, i, j
  real mini,maxi
  integer  :: t1, t2, rate
  real, allocatable :: tmp (:,:,:)
  real, allocatable :: work(:)
!  double precision :: mean
  real :: mean
  real :: slange

  call system_clock(count_rate=rate)
  call system_clock(t1)
  gridsize=512
  unit=40

  allocate( tmp(gridsize, gridsize, gridsize), work(gridsize))

  open(unit=unit,file='T.out',status='old',access='stream',&
       form='unformatted',action='read')
  read(unit=unit) tmp

  close(unit=unit)

  mini = minval(tmp)
  maxi = maxval(tmp)

!  mean = sum(tmp)/gridsize**3
!  mean = sum(real(tmp, kind=kind(1.d0)))/real(gridsize**3, kind=kind(1.d0))
  mean = 0.d0
  do j=1,gridsize
    do i=1,gridsize
      mean = mean + slange('1', gridsize, 1, tmp(:,i,j),gridsize, work)
    enddo !i
  enddo !j
  mean = mean / gridsize**3

  print *, mini, maxi, mean
  call system_clock(t2)
  print *,real(t2-t1)/real(rate)

end program

这在矩阵列上使用单精度矩阵 1-norm SLANGE。运行时间甚至比使用单精度数组函数的方法更快 - 并且没有显示精度问题。

【讨论】：

为什么混合输入和计算会减慢速度？他们都必须阅读整个文件，这将是瓶颈。如果操作系统确实预读，则 Fortran 代码不必为 I/O 等待太多时间。

@Barmar 您仍然需要每次都检查数据是否在缓存中的函数调用开销和逻辑。

【参考方案2】：

numpy 更快，因为您在 python 中编写了更高效的代码（并且大部分 numpy 后端是用优化的 Fortran 和 C 编写的）并且在 Fortran 中编写的代码效率非常低。

查看您的 python 代码。您一次加载整个数组，然后调用可以对数组进行操作的函数。

查看您的 fortran 代码。您一次读取一个值并使用它执行一些分支逻辑。

您的大部分差异是您在 Fortran 中编写的碎片化 IO。

您可以像编写 python 一样编写 Fortran，您会发现它运行得更快。

program test
  implicit none
  integer :: gridsize, unit
  real :: mini, maxi, mean
  real, allocatable :: array(:,:,:)

  gridsize=512
  allocate(array(gridsize,gridsize,gridsize))
  unit=40
  open(unit=unit, file='T.out', status='old', access='stream',&
       form='unformatted', action='read')
  read(unit) array    
  maxi = maxval(array)
  mini = minval(array)
  mean = sum(array)/size(array)
  close(unit)
end program test

【讨论】：

以这种方式计算的平均值是否与numpy的.mean调用获得相同的精度？我对此有些怀疑。

@Bakuriu 不，它没有。请参阅 Alexander Vogt 的回答和我对该问题的编辑。