我从一个模拟(用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例程中做了一些非常愚蠢的事情,让它花了这么长的时间?
编辑:
要回答评论中的问题:
iso_fortran_env
。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
如何避免这种情况呢?
感谢你到目前为止所做的一切。
https://stackoverflow.com/questions/33723771
复制相似问题