numpyは私のFortranルーチンよりもずっと高速でしょうか？

Question

シミュレーション（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) 

は今、私は、これは当然の速くなると予想したが、私は本当に吹き飛ばされました。同じ条件下では1秒もかかりません。私のFortranルーチンが見つけたものからの平均（128ビット浮動小数点でも動いていたので、何とかそれをもっと信頼しています）が、7番目の有意な桁にすぎません。

どのようにnumpyを高速にすることができますか？これらの値を見つけるために配列のすべてのエントリを調べなければならないのですよね？私はFortranのルーチンで非常に愚かなことをしていますか？

EDIT：

コメントでの質問に答えるために：

• はい、また私は、32ビットおよび64ビットの浮動小数点数とFortranのルーチンを実行しましたが、それは、パフォーマンスに影響を及ぼしませんでした。
• 私は128ビット浮動小数点を提供するiso_fortran_envを使用しました。
• 32ビット浮動小数点を使用すると、私の平均値はかなり小さくなります。したがって、精度は本当に問題です。
• 私は両方のルーチンを別々のファイルで異なる順序で実行しました。だから、キャッシュは私が推測している比較の中で公正であったはずですか？
• 私は実際に開いたMPを試みましたが、同時に異なる位置でファイルから読み込みました。あなたのコメントと答えを読んだら、これは本当にばかげた音となり、日常生活にもかなり時間がかかりました。私はそれを配列操作で試してみるかもしれませんが、多分それは必要ではないでしょう。
• ファイルは実際にはサイズが1/2Gです。これはタイプミスです。ありがとうございました。
• ここで配列の実装を試みます。

EDIT 2：

私はその答えで提案されているものを@Alexanderフォークトと@casey実装し、それが早くnumpyようですが、@Luaanが、私は可能性があります指摘したように、今私は、精度の問題を抱えています取得する。 32ビット浮動小数点配列を使用すると、sumで計算された平均は20％オフです。 Doing

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

問題を解決しますが、計算時間が増加します（ただし、それほど顕著ではありません）。 この問題を回避するより良い方法はありますか？私は、ファイルから直接ダブルスにシングルを読む方法を見つけることができませんでした。 そして、どうすればnumpyはこれを避けますか？

これまでのお役に立てていただきありがとうございます。

Answer 1

あなたのFortranの実装は二つの主要な欠点を被る：

• あなたはIOと計算をミックス（およびエントリによるファイルのエントリから読み取ります）。
• ベクトル/行列演算は使用しません。

これ実装はあなたと同じ操作を実行して、私のマシン上で速く20倍である：

アイデアは一度に1列 tmpにファイル全体を読み込むことです

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 

。次に、配列の関数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を単精度行列を使用します。実行時は、単精度配列関数を使用するアプローチよりも高速であり、精度の問題は示されません。あなたはpythonではるかに効率的なコードを書いた（とnumpyのバックエンドの多くは最適化されたFortranとCで書かれている）およびFortranでひどく非効率的なコードので

Answer 2

numpyのは速いです。

は、あなたのPythonコードを見てください。アレイ全体を一度にロードしてから、配列で操作できる関数を呼び出します。

はあなたのFortranコードを見てください。あなたは一度に1つの値を読んで、それと何らかの分岐ロジックを行います。

あなたの不一致の大半は、あなたが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 test