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怎么能这么快？我的意思是，你必须查看数组的每个条目才能找到这些值，对吧？我是否在我的Fortran例程中做了一些非常愚蠢的事情，让它花了这么长的时间？
编辑：
要回答评论中的问题：
- 是的，我用32位和64位浮点数运行了Fortran例程，但它对性能没有影响。
- 我使用了提供128位浮点数的
iso_fortran_env。 - 使用32位浮点数我的意思是，我的意思是不太准确，所以精度确实是一个问题。
- 我以不同的顺序在不同的文件上运行这两个例程，所以缓存应该是公平的，我猜？<代码>H213<代码>H114我实际上尝试打开MP，但同时在不同的位置读取文件。在阅读了你的评论和答案后，这听起来真的很愚蠢，而且它也让例程花费了更长的时间。我可能会在数组操作上尝试一下，但也许那不会是necessary.
- The文件的大小实际上是1/2G，这是一个拼写错误，谢谢。
- 我现在将尝试数组实现。
EDIT 2:
我实现了@Alexander Vogt和@casey在他们的答案中建议的，它和numpy一样快，但现在我有一个精度问题，@Luaan指出我可能会得到。使用32位浮点数组，sum计算的平均值为20%。正在做什么
...
real,allocatable :: tmp (:,:,:)
double precision,allocatable :: tmp2(:,:,:)
...
tmp2=tmp
mean=sum(tmp2)/size(tmp)
...解决了这个问题，但增加了计算时间(不是很多，但很明显)。有没有更好的方法来解决这个问题？我找不到一种方法直接从文件中读取单元格到双元格。numpy如何避免这种情况呢？
感谢你到目前为止所做的一切。
回答 2
Stack Overflow用户
发布于 2015-11-16 04:07:34
您的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范数SLANGE。运行时甚至比使用单精度数组函数的方法更快-并且没有显示精度问题。
Stack Overflow用户
发布于 2015-11-16 04:18:32
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 testhttps://stackoverflow.com/questions/33723771
复制相似问题