fortranにおいてpythonをメインプログラムにするには？

前提: fortranの仕事はmodule内にデータを貯めること

fortranは数値計算に適しているが、pythonのもつ自由度がない。将来的にはpythonを親に据えたコードに書き換えていく必要がある。上位でのMPI制御はpython側で行うべき。

fortranで巨大プログラムを書くのはやめたほうがいいーmoduleによるsingletonしかできないものだと悟るべき。すなわち、fortranの仕事を「fortranのmodule内にデータを貯める」という操作だけにとどめるべきである.fortranのmodule内subroutineは指示のみを引数として受け取り、その指示に従ってデータをmodule内に貯めこむという操作のみを行う(protectedで書く)。データは他のmoduleやファイルから読み込む。これを繰り返して計算を行うことになる。（将来的にはfortran code内でMPI関数のコールはしないほうがいいと思われる。）

fortranをsingletonパターンを要請する「上記のことしかできない数値計算専用プログラミング言語」として見るのならそんなに悪くはない。むしろ特化している分、ムダなお祈り文を書かなくてすむ。正直、fortranはその路線で進化すべきであり、汎用言語として進化する必要はないと思う。

pythonからfortranを呼び出すサンプル

こういうことになると、およそ単純なsubroutine callの繰り返しで計算が進行していくことになる。ただしpythonからライブラリ化したfortranのsubroutineを呼ばないといけない。またMPIはpython側で起動しcommunicatorをfortranに渡さないといけない。

そのサンプルを書いてみた.メインプログラムhello.pyは非常に簡単になる。

# main program hello.py
from setcomm import *
ecaljF = setcomm("ecaljF.so")
callF( ecaljF.hello )
callF( ecaljF.hello )
callF( ecaljF.hello3, [True, 34, 8.2, 9.3] )

このサンプルでは,fortranのsubroutine helloとsubroutine hello3を呼び出している。以下がそのソースコードである。これらのsubroutineはインプット引数しか許さないし、pythonへの返り値はなし。fortranコードの書換として必要なのは「fortran側ではbind(C)をつけること。MPIのsize,rank,communicatorはmodule m_commの変数を読む」という点のみであるので単純である。関数呼び出しのとき、module内関数はグローバルにアクセスできるものとなっている点には注意がいる（useみたいなものを指定しない）。

#fortran subroutines 
module m_test1
  use m_comm,only: size,rank
  !include "mpif.h"
  integer,protected:: iii,n
  logical:: init=.true.
  real(8):: a=3.21
  complex(8),allocatable:: ccc(:,:)
  private
  !public hello
contains
  subroutine hello() bind(C)
    implicit none
    complex(8):: cc
    if(init) then
       print *,'rank/size=',rank,size, "Hello World111init"
       init=.false.
       a=1.23d0
       n=1000
       cc=rank
       allocate(ccc(n,n),source=cc)
    else
      print *,'rank/size=',rank,size, "Hello World111",a,sum(ccc)/n/n
    endif
  end subroutine hello
end module m_test1
module m_test2
contains
 subroutine hello3(lll,intx,r1,r2) bind(C)
   use m_comm,only: comm,size,rank
   implicit none
   include "mpif.h"
   integer :: ierr,intx
   logical:: lll
   real(8):: r1,r2
   call MPI_barrier(comm,ierr)
   if(mod(rank,2)==0) print *,'rank/size=',rank,size, "Hello World333 even", lll,intx,'sum=',r1+r2
   if(mod(rank,2)==1) print *,'rank/size=',rank,size, "Hello World333 odd ", lll,intx,'sum=',r1+r2
 end subroutine hello3
end module m_test2

hello()ではテストのためallocationなどmodule変数へのデータ格納を試している。

上の、python-fortranの結合のからくりにおいては関数setcomm("ecaljF.so")がミソになる。これは以下の2つのファイルsetcomm.pyとm_comm.pyで定義されている。

# setcomm.py
from ctypes import *
import numpy as np
from mpi4py import MPI

def setcomm(aaa): #Pass communicator to module m_comm.f90
    comm = MPI.COMM_WORLD
    comm = comm.py2f()
    eee = np.ctypeslib.load_library(aaa,".")
    eee.setcomm.argtypes = [ POINTER(c_int32) ]
    #eee.setcomm.restype  = c_void_p
    eee.setcomm(c_int32(comm))
    return eee

def callF(eee,mmm=[]):
    argtypess=[]
    data=[]
    for ii in  mmm:
        #print(ii,type(ii))
        if(type(ii)==type(1)): 
            argtypess.append( POINTER(c_int32) )
            data.append(c_int32(ii))
        elif(type(ii)==type(1.0)):
            argtypess.append( POINTER(c_double) )
            data.append(c_double(ii))
        elif(type(ii)==type(True)):
            argtypess.append( POINTER(c_bool) )
            data.append(c_bool(ii))
    eee.argtypes= argtypess
    if(len(argtypess)==0):
        eee()
    else:
        eee(*data)
    return

# m_comm.f90
module m_comm
  integer,protected:: comm,size,rank
contains
  subroutine setcomm(commin) bind(C)
    implicit none
    include "mpif.h"
    integer :: ierr,n,commin
    call MPI_Comm_size(commin, size, ierr)
    call MPI_Comm_rank(commin, rank, ierr)
    comm=commin
  end subroutine setcomm
end module m_comm

汎用コードであるので、fortranユーザーはいちいちsetcomm.pyとm_comm.pyをいじる必要はない。function呼び出しなどにも拡張できる。正直、十分に単純なのでコンパイラー＆リンカーがトラブルを起こすことは少ないと思う。

一回目のコールのcallF( ecaljF.hello )でモジュール変数を貯めており
２回目のcallF( ecaljF.hello )のコールでそれを参照して違う結果を与えることを考えることができる。
hello.pyはfortranでも簡単に書き換えれて

#hello.f90
program hello
use m_comm,only: setcomm
implicit none
include "mpif.h"
integer:: info
call MPI_Init( info )
call setcomm(MPI_COMM_WORLD)
call hello()
call hello()
call hello3(True, 34, 8.2d0, 9.3d0)
end

となる。fortranで完結させて動作チェックした上でpython化しないといけないので、この単純さとhello.pyとの間の単純な対応関係は重要である。
このレベルであれば、hello.f90をhello.pyへ自動書き換えを行うスクリプトも単純に書ける思われる。

ソース

https://github.com/tkotani/PythonFortran/blob/master/hello.py
に置いた。READMEがあるが、

 mpif90 -shared -fPIC -o ecaljF.so *.f90 #ecaljF.soを作る
 mpiexec -n 4 python3 ./hello.py         #4coreでMPI実行　　　　　　

で実行できる。

そもそもpythonがfortranに指示を送ってfortranのモジュール内にデータを貯め込む、という設計思想を守るべきなので、python-fortranの間で複雑なデータのやり取りはさせない。pythonにプロットさせたければ、fortranでファイルを書かせてそれを受け取って書かせるので良い。そんなところでのバイナリレベルでの結合に時間を費やすのは時間のムダだしトラブルのもと。とくにfortran側を書き換えたくない。

ただそうは言っても教条主義的なのが一番いかんことだし、もう少し改良してfortranのreal(8),integerの一次元配列ぐらいは同じやり方で渡せるようにしたい。さらには

# main program hello.py
from setcomm import *
ecaljF = setcomm("ecaljF.so")
callF( ecaljF.hello )
callF( ecaljF.hello )
callF( ecaljF.hello3, [True, 34, 8.2, 9.3, array1input],[out1,out2] )

などと書けるとうれしいな。どうやればいいのかな？(１．fortranのほうはcomm参照とbind(C)以外書き換えたくない。２．pythonと結合しないと動かないようなfortranコードでは困るし、そのために大きく書き換えないといけないのでは困る。３．返り値out1,out2も受け取れるならその方が良い。）

いいのができたら教えてほしいです