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中。 然后,我可以直接使用数组上的函数MAXVALMINVALSUM


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

 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