numpy怎么能比我的Fortran程序快得多呢？
我从一个模拟（用Fortran编写）得到一个表示温度分布的512 ^ 3数组。 该数组存储在大小约为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怎么能这么快？ 我的意思是你必须看看数组的每一个条目来find这些值，对吧？ 我在Fortran程序中做了一些非常愚蠢的事情，因为它需要更长的时间？
编辑：
要回答评论中的问题：
- 是的，我也使用32位和64位浮点运行Fortran程序，但对性能没有影响。
- 我使用了提供128位浮点的
iso_fortran_env。 - 使用32位浮点数我的意思是closures不less，所以精确度是一个问题。
- 我以不同的顺序在不同的文件上运行了两个例程，所以caching应该是公平的比较，我猜？
- 我其实试过打开MP，但是同时从不同位置的文件中读取。 看完你的评论和回答，这听起来真的很愚蠢，这也使得例行公事也变得更长。 我可能会尝试一下arrays操作，但也许这不是必须的。
- 这些文件实际上是1 / 2G，这是一个错字，谢谢。
- 我将尝试现在的数组实现。
编辑2：
我在他们的答案中实现了@Alexander Vogt和@casey的build议，它的速度和numpy一样快，但是现在我有一个精确的问题，就像@Luaan指出的那样。 使用32位浮点数组， sum计算的平均值为20％。 干
... real,allocatable :: tmp (:,:,:) double precision,allocatable :: tmp2(:,:,:) ... tmp2=tmp mean=sum(tmp2)/size(tmp) ...
解决这个问题，但增加了计算时间（不是很多，但明显）。 有没有更好的方法来解决这个问题？ 我找不到直接从文件中读取单曲的方法来双打。 而numpy如何避免这种情况？
感谢所有迄今为止的帮助。
您的Fortran实现存在两个主要缺陷：
- 混合IO和计算（并通过条目从文件条目中读取）。
- 你不使用vector/matrix操作。
这个实现和你一样执行相同的操作，在我的机器上速度提高了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 ，元素相加比arrays操作要好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
这在matrix列上使用单精度matrix1-范数SLANGE 。 运行时比使用单精度数组函数的方法更快 – 并且不显示精度问题。
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
