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

Question

我得到一个 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 又是如何避免这种情况的呢？

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

Answer 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。运行时间甚至比使用单精度数组函数的方法更快 - 并且没有显示精度问题。

Answer 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