## 为什么Numpy比我的Fortran程序快得多？内容来源于 Stack Overflow，并遵循CC BY-SA 3.0许可协议进行翻译与使用

• 回答 (2)
• 关注 (0)
• 查看 (22)

```  integer gridsize,unit,j
real mini,maxi
double precision mean

gridsize=512
unit=40
open(unit=unit,file='T.out',status='old',access='stream',&
mini=tmp
maxi=tmp
mean=tmp
do j=2,gridsize**3
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)```

```    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)```

```...
real,allocatable :: tmp (:,:,:)
double precision,allocatable :: tmp2(:,:,:)
...
tmp2=tmp
mean=sum(tmp2)/size(tmp)
...```

### 2 个回答

• 您可以混合IO和计算(并按条目从文件条目读取)。
• 您不使用向量/矩阵操作。

```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',&

close(unit=unit)

mini = minval(tmp)
maxi = maxval(tmp)
mean = sum(tmp)/gridsize**3
print *, mini, maxi, mean

end program```

`mean = sum(real(tmp, kind=kind(1.d0)))/real(gridsize**3, kind=kind(1.d0))`

`-O3`，元素级加法比数组运算执行约3%的运算。下面是使用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',&

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```

```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',&