fortranをもっと便利に使いたい

fortranに関して、2019年時点ではf2008規格へのコンパイラ対応も充実し、オブジェクト指向も記述できるようになっている。f77の書き方はオブジェクト指向を含む処理より計算速度的には有利なので、両者をうまく使えば、ScipyやOctave並みの使い易さを持ちつつ、計算速度の非常に高速なライブラリができるのではないか。f77の資産¹を活用する場面でひと工夫したので、具体的な事例を記録しておく。

具体事例

ScipyやOctaveの説明書を読んでいると、fortranのライブラリを参照する説明がしばしば登場する。fortranユーザーなら、それらのライブラリを直接使えば事足りるとはいえ、それらはf77の様式で書かれていることが多い。そこで、これを呼び出すオブジェクトを作り、不必要な引数は省略可能にしつつ、必要最小限の引数で便利な関数を使えるようにしたい。そうすれば、ScipyやOctaveと同じような感覚で数値計算ライブラリを呼び出すことができる。以下では、例として非線形方程式の解法に着目する。

Octaveの例

例えば、Octaveでは非線形方程式

\boldsymbol{f}(\boldsymbol{x})=\boldsymbol{0}

fsolve(FCN,X0,OPTIONS)

ここに、FCNはユーザ定義関数、X0はXの初期値、OPTIONSは計算のパラメータ等を指定することができるが、デフォルト値で良い場合はOPTIONSを省略可能である。すなわち、ベクトルxを与えた時の関数ベクトルfさえ計算できればOctaveで何らかの解を得ることができる。Octaveのマニュアルによると、fsolve関数はfortranのMINPACKに含まれるhybrdに基づくらしいので、次にこれを見てみる。

f77の例

MINPACKのhybrd関数を見てみると、Octaveよりはるかに多くの引数が必要である。

   subroutine hybrd(fcn,n,x,fvec,xtol,maxfev,ml,mu,epsfcn,diag,&
                    mode,factor,nprint,info,nfev,fjac,ldfjac,r,lr,&
                    qtf,wa1,wa2,wa3,wa4)
      integer n,maxfev,ml,mu,mode,nprint,info,nfev,ldfjac,lr
      double precision xtol,epsfcn,factor
      double precision x(n),fvec(n),diag(n),fjac(ldfjac,n),r(lr),&
                       qtf(n),wa1(n),wa2(n),wa3(n),wa4(n)
      external fcn

さすがに引数が多いのか、MINPACK自体にも、引数の少ないバージョンが用意されていて、次のようなインターフェイスとなっている。

    subroutine hybrd1(fcn,n,x,fvec,tol,info,wa,lwa)
      integer n,info,lwa
      double precision tol
      double precision x(n),fvec(n),wa(lwa)
      external fcn

hybrdよりもhybrd1はかなり引数の数が少なくなっているが、それでもOctaveよりも多く、f77なので省略可能な引数はない。Octaveでは、おそらく最初の4つの引数のみをユーザに指定させ、その他はデフォルト値としてOPTIONSに含めているものと思われる。
Octaveなどの言語を知ってしまった身としては、できればtol以降の引数はユーザが細かく制御可能な状態のまま、通常はコンピュータに自動で宜しきにやってもらいたい。そこで、fortran2003以降に導入されたオブジェクト指向を使ってこれを実現する。

オブジェクト指向を活用して使いやすく

fortran2003からは、データとメソッドを持つ構造体を定義できて、さらにそれを継承する拡張型を定義できる。さらに、オブジェクト内部の公開属性を指定できるので、カプセル化も可能である。
具体的には、Modern Fortran Explainedという本がとても詳しい。日本語ではFortranからみたオブジェクト指向などが分かり易いと思う。

fortranで非線形方程式を解くオブジェクトを作る

hybrd1をf77でコールする時に必要な情報を含む、次のような構造体を定義する。

 type,public :: hybrdtype 
    integer :: info=-1 !<(out) hybrd1からの情報
    integer :: n    !<(in)関数値fと変数xの要素数
    integer :: lwa  !<hybrd1に渡されるwork arrayの要素数
    double precision :: tol=1d-10 !<(in) 誤差制御のデフォルト値
    double precision,dimension(:),allocatable :: x !<(inout)変数
      !<インプットは初期推定値、アウトプットはhybrd1による最終推定値。
    double precision,dimension(:),allocatable :: fvec !<(out)関数値
    double precision,dimension(:),allocatable :: wa !<work array
    procedure(),nopass,pointer::fcn=>null() !<fvecを計算する関数ポインタ
    contains
      final :: final_hybrd !< デストラクタ
      procedure :: setup => setup_hybr_common !<ユーザーサプライ関数や変数のセット
      procedure :: run   => run_hybrd   !<minpackのhybrd1関数を実行
  end type

contains以下のsetup手続きは、構造体内部の配列やポインタを定義するためのもので、次のインターフェイスとした。

subroutine setup_hybr_common(self,fcn,n,tol,xinit)
  class(hybrdtype),intent(inout) :: self
  integer,intent(in) :: n
  procedure() :: fcn 
  double precision,optional :: tol
  double precision,dimension(n),optional :: xinit

なお、run_hybrdは、定義したオブジェクトを用いてf77のhybrd1ルーチンをコールしているだけである。このようなオブジェクトを定義しておけば、次のような簡潔なコードで済むため、利用するときに型を宣言してセットアップして実行するだけとなり、引数の多さに悩まされない。

program test
  use hybrdmodule !hybrdtypeの定義を記述したモジュール
  implicit none
  type(hybrdtype) :: hybrd
  integer :: n=2
  !関数の定義と方程式&未知数の数、さらに初期値を設定
  call hybrd%setup(extfcn,n,xinit=[0d0,0d0]) 
  call hybrd%run() !hybrd1をcall
  print*,hybrd%x(:) !解を表示
end program

ちなみに、extfcnは例えば次のような関数を定義すれば良い。

  subroutine extfcn(n,x,fvec,iflag)
    integer :: n,iflag
    double precision::x(n),fvec(n)
    fvec(1)=-2d0*x(1)*x(1)+3d0*x(1)*x(2)+4d0*sin(x(2))-6d0
    fvec(2)=3d0*x(1)*x(1)-2d0*x(1)*x(2)**2d0+3d0*cos(x(1))+4d0
  end subroutine

おわりに

今回は、MINPACKを例にf77のコードを使いやすくする例を記した。この方法はODEPACKなどにも適用できるため、よく使うルーチンは同様の方法でモジュールにまとめている。なお、f77のルーチン毎にfortranの構造体を定義するのではなく、似たルーチンはオブジェクトの拡張やメソッドのオーバーライドを活用することもできる。
さらに、色々とモジュールを作っていると、ひとまとめにしてパソコンにインストールしたくなり、configure && make && make install でライブラリとして使えるようにしている。fortranの情報は少ないので、いつか時間があればfortranのconfigureスクリプトの作り方をメモしたい。

そもそも計算アルゴリズムを近年のコンピュータ向けに書き直すという考えは脇に置いておく。 ↩
xもfも,要素が2つ以上あってよい。
を解く時に次の関数を使うことができる²。 ↩