読者です 読者をやめる 読者になる 読者になる

趣味で計算流砂水理

Computational Sediment Hydraulics for Fun

pythonのfortranによる高速化

R

※ifortで通ったので修正しました。

I/Oはpythonで書くけど、科学技術計算の部分はfortranに任して高速化してしまおうという話です。以下のページをほぼそのまま使用しました。

PythonからFortranのサブルーチンを呼ぶ。 - いぐにさんにっき

windows 環境で fortran と C++ を連携させる(intel fortranの場合) | Moonmile Solutions Blog

ちなみにnumbaによる高速化も試しましたが、今回のプログラムでは全然早くなりませんでした。どうもリスト内包標記が駄目っぽいです。それだとpythonの良さがなくなってしまうので使えないですね。

速度比較

  • モデルは例のsectionQ2Dモデルです。
  • 収束計算はガウス・ザイデルです。(これをADIとかに変えれば良いという話もありますが。)

  • とりあえず速度の比較です。

    • python:231.958 sec
    • gfortran:1.222 sec
    • ifort:0.706 sec

衝撃の速さです。なお、イテレーションの回数は同じです。

つまずいた点

  • 計算条件等はpythonで設定して、配列のままfortranに渡すことを試みました。環境依存は結構ありそうです。今回は以下の環境です。

    • windows7 64bit
    • python3.5
    • gfortran 4.8.1
    • ifort 2013
  • つまずいた点は以下です。

    • gfortranの場合、pythonから呼ぶ際にsubroutin名の末尾にアンダースコアが必要。
    • ifortの場合、pythonから呼ぶ際にsubroutin名がすべて大文字。(32bitの場合はアンダースコアが必要らしい。環境がないので試してません。)
    • 多次元配列が逆になる。←CとFortranが逆のことを考えると当たり前ですが、普通に通っていたので全然気付きませんでした。 - gfortranの場合、functionは使えない。←subroutineだけっぽいです。
    • Fortranの動的配列は使えない。←配列の大きさを別途渡す必要があります。

ソースコード

python

def sectionq2d3():
    un = np.array([[ini_velocity for k in range(nz)] for j in range(ny)])
    unnew = np.zeros((ny,nz))

    def NumericalSimulation():
        for j in range(0,ny):
            for k in range(0,nz):
                if j == 0 : # Wall Boundary
                    apy = (nyut2[k]+nyut2[k])*dz**2
                    ayp = dz**2*nyut2[k]*un[j+1,k]
                    aym = 0
                elif j == ny-1 : # Wall Boundary
                    apy = (nyut2[k]+nyut2[k])*dz**2
                    ayp = 0
                    aym = dz**2*nyut2[k]*un[j-1,k]
                else: # Normal
                    apy = (nyut2[k]+nyut2[k])*dz**2
                    ayp = dz**2*nyut2[k]*un[j+1,k]
                    aym = dz**2*nyut2[k]*un[j-1,k]

                if k == 0 : # Bed Boundary
                    apz = (0.5*nyut2[k+1]+nyut2[k]+0.5*nyut2[k])*dy**2
                    azp = dy**2*0.5*(nyut2[k]+nyut2[k+1])*un[j,k+1]
                    azm = 0
                elif k == nz-1 : # WaterSurface Boundary: Slip Condition
                    apz = 0.5*(nyut2[k]+nyut2[k-1])*dy**2
                    azp = 0
                    azm = dy**2*0.5*(nyut2[k]+nyut2[k-1])*un[j,k-1]
                else: # Normal
                    apz = (0.5*nyut2[k+1]+nyut2[k]+0.5*nyut2[k-1])*dy**2
                    azp = dy**2*0.5*(nyut2[k]+nyut2[k+1])*un[j,k+1]
                    azm = dy**2*0.5*(nyut2[k]+nyut2[k-1])*un[j,k-1]

                un[j,k] = (9.8*ib*dy**2*dz**2+ayp+aym+azp+azm)/(apy+apz)

    for num in range(1,10000):
        NumericalSimulation()
        if np.max(np.abs(unnew-un)) > 0.000001 : unnew[:,:] = un
        else :
            print(num)
            break

    return un

fortran

python側でdllを呼ぶ。

from ctypes import *

def sectionq2d4(fortType):
    if fortType == 'gfort' :
        f = np.ctypeslib.load_library("sq2dFortranGfort.dll",".")
        f.sectionq2d4_.argtypes = [
                np.ctypeslib.ndpointer(dtype=np.float64),
                np.ctypeslib.ndpointer(dtype=np.float64),
                POINTER(c_int64),POINTER(c_int64),
                POINTER(c_double),POINTER(c_double),POINTER(c_double)
                ]
        f.sectionq2d4_.restype = c_void_p
    elif fortType == 'ifort' :
        f = np.ctypeslib.load_library("sq2dFortranIfort.dll",".")
        f.SECTIONQ2D4.argtypes = [
                np.ctypeslib.ndpointer(dtype=np.float64),
                np.ctypeslib.ndpointer(dtype=np.float64),
                POINTER(c_int64),POINTER(c_int64),
                POINTER(c_double),POINTER(c_double),POINTER(c_double)
                ]
        f.SECTIONQ2D4.restype = c_void_p

    un = np.array([[ini_velocity for k in range(nz)] for j in range(ny)],dtype=np.float64)
    nyut21 = np.array(nyut2,dtype=np.float64)
    ny2 = byref(c_int64(ny))
    nz2 = byref(c_int64(nz))
    dy2 = byref(c_double(dy))
    dz2 = byref(c_double(dz))
    ib2 = byref(c_double(ib))

    if fortType == 'gfort' : f.sectionq2d4_(un,nyut21,ny2,nz2,dy2,dz2,ib2)
    elif fortType == 'ifort' : f.SECTIONQ2D4(un,nyut21,ny2,nz2,dy2,dz2,ib2)

    return un

fortranのコードです。

subroutine sectionq2d4(un, nyut2, ny, nz, dy, dz, ib)
! ifortの場合、以下の1文が必要。gfortranの場合は不要。
    !DEC$ ATTRIBUTES DLLEXPORT :: sectionq2d4
    implicit none
    integer(8), INTENT(IN) :: ny, nz
    real(8), INTENT(IN) :: dy, dz, ib
    real(8), INTENT(IN) :: nyut2(0: nz-1)
    real(8), INTENT(INOUT) :: un(0: nz-1, 0: ny-1)
    integer(8) :: num
    real(8) :: unnew(0: nz-1, 0: ny-1)
    real(8) :: tmp(0: nz-1, 0: ny-1)
    integer(8) :: k, j
    real(8) :: apy,ayp,aym,apz,azp,azm

    unnew = 0.0

    do num = 1, 1000000
        call NumericalSimulation

        if (maxval(abs(unnew-un)) > 0.000001) then
            unnew = un
        else
            print *,num
            exit
        end if
    end do
CONTAINS
    subroutine NumericalSimulation()
        do j = 0, ny-1
            do k = 0, nz-1
                if (j == 0) then  !Wall Boundary
                    apy = (nyut2(k)+nyut2(k))*dz**2
                    ayp = dz**2*nyut2(k)*un(k,j+1)
                    aym = 0
                else if (j == ny-1) then ! Wall Boundary
                    apy = (nyut2(k)+nyut2(k))*dz**2
                    ayp = 0
                    aym = dz**2*nyut2(k)*un(k,j-1)
                else  ! Normal
                    apy = (nyut2(k)+nyut2(k))*dz**2
                    ayp = dz**2*nyut2(k)*un(k,j+1)
                    aym = dz**2*nyut2(k)*un(k,j-1)
                end if
                if (k == 0) then ! Bed Boundary
                    apz = (0.5*nyut2(k+1)+nyut2(k)+0.5*nyut2(k))*dy**2
                    azp = dy**2*0.5*(nyut2(k)+nyut2(k+1))*un(k+1,j)
                    azm = 0
                else if (k == nz-1) then ! WaterSurface Boundary: Slip Condition
                    apz = 0.5*(nyut2(k)+nyut2(k-1))*dy**2
                    azp = 0
                    azm = dy**2*0.5*(nyut2(k)+nyut2(k-1))*un(k-1,j)
                else ! Normal
                    apz = (0.5*nyut2(k+1)+nyut2(k)+0.5*nyut2(k-1))*dy**2
                    azp = dy**2*0.5*(nyut2(k)+nyut2(k+1))*un(k+1,j)
                    azm = dy**2*0.5*(nyut2(k)+nyut2(k-1))*un(k-1,j)
                end if
                un(k,j) = (9.8*ib*dy**2*dz**2+ayp+aym+azp+azm)/(apy+apz)
            end do
        end do
    end subroutine
end subroutine

fortranコンパイルです。

gfortran -shared -o sq2dFortran.dll sq2dFortran.f90
ifort /nologo /dll /iface:cvf -o sq2dFortranIfort.dll sq2dFortran.f90

相当使えそうです。やっぱりfortranはサルでも速いコードがかけるので便利ですね。