numpy如何比我的Fortran例程快得多?-问答-阿里云开发者社区-阿里云

开发者社区> 问答> 正文

numpy如何比我的Fortran例程快得多?

保持可爱mmm 2020-02-08 14:35:45 71

我从仿真中得到了一个代表温度分布的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,gridsize3 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/gridsize3 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例程中做了一些非常愚蠢的事情,以使其花费了更长的时间? 问题来源于stack overflow

存储 Python
分享到
取消 提交回答
全部回答(1)
  • 保持可爱mmm
    2020-02-08 14:36:01

    您的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)/gridsize3 ! mean = sum(real(tmp, kind=kind(1.d0)))/real(gridsize3, 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 这SLANGE在矩阵列上使用单精度矩阵1-范数。运行时间甚至比使用单精度数组函数的方法快-并且没有显示精度问题。

    0 0
开发与运维
使用钉钉扫一扫加入圈子
+ 订阅

集结各类场景实战经验,助你开发运维畅行无忧

推荐文章
相似问题
推荐课程