为什么Numpy比我的Fortran程序快得多?

内容来源于 Stack Overflow,并遵循CC BY-SA 3.0许可协议进行翻译与使用

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

我得到了一个512^3数组,它代表了一个模拟温度分布(用Fortran编写)。数组存储在一个二进制文件中,该文件大小约为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位左右。

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

解决了这个问题,但增加了计算时间(不是很多,而是明显的)。有没有更好的方法来解决这个问题?如何numpy避开这个?

提问于
用户回答回答于

您的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一步一步。然后,我可以使用这些函数MAXVALMINVAL,和SUM直接在数组上。

对于精度问题:只需使用双精度值,并执行动态转换,如

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',&
       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在矩阵列上。运行时甚至比使用单精度数组函数的方法还要快,并且没有显示精度问题。

用户回答回答于

您可以使用与编写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

扫码关注云+社区