R で Fortran のサブルーチンを使う：R と Fortran の結合方法

R で Fortran のサブルーチンを使う：R と Fortran の結合方法
参照:
http://d.hatena.ne.jp/rontas/20110621/1308614194
0. gfortran をインストールする
http://gcc.gnu.org/wiki/GFortranBinaries の上の，例えば，
「Download the latest installer
(dated 2012-12-10)」をクリックすると，gfortran-windows-20121210.exe をダウンロード
されるので，それをインストールする． MinGW for Win64 の gfortran については，未確
認．gfortran の使い方，実行例は，
http://www.esst.kyushu-u.ac.jp/~space/kumashiro/fortran.html
などを参照．
1. Fortran サブルーチンを用意する
以下のような Fortran 形式のサブルーチンを，秀丸やメモ帳などで test.for という名前を付
けて保存する．このとき，拡張子を f90 にすると fortran90 のファイル形式となる(for ファ
イルと異なり，最初のタブがいらない)．
subroutine test(a, b)
double precision a, b
b = 1000.0*a + 100.0*a + 10.0*a + a
end
3. gfortran で Fortran サブルーチンをコンパイルする
コマンドプロンプトを立ち上げ，（次ページの図のように）「cd ファイル場所」を使い，
Fortran ファイルが存在するフォルダに作業フォルダを移動したうえで，以下のような命令
を打ち込む：
gfortran --shared -o test.dll test.for
※ここで、重要なのは、サブルーチンの名前と保存する Fortran コードの名前、作成する
dll の名前を一致させる必要があるということである。また、このコードを保存するディレ
クトリの名前は全て小文字にした方がうまくいく場合が多い。またコンパイルは、
R CMD SHLIB test.f90 などとしてもよいが、R 側の設定が少々必要である。
4. R にロードする
ここからは、R コンソール上での作業になる。以下
> dyn.load("test.dll")
のように打ち込めば、R 上に test.dll がロードされる。きちんとロードされているかどうか
確かめるには以下のようにする：
> is.loaded("test") または getLoadedDLLs()
5. R から呼び出す関数を作る
4 に続いて，R 上で，以下のように用いれば
> a=1
> .Fortran("test",arg1=as.double(a),arg2=double(1))
$arg1
[1] 1
$arg2
[1] 1111
という結果が返され，$arg2 の結果 b (a=1 のときの)が返され、Fortran のサブルーチンを
使うことができるようになる。
6. その他：起こりえるエラーと対処
その他に，dll のロード時や予想外のエラーが出れば，以下を確認すること

gfortran を 32-bit でインスールすれば，R も 32-bit バージョンで使う必要があると思
う．

Window 7 より前の PC では「libgcc_s_dw2-1.dll」も予め（このファイルはネット検
索すれば，どこかに落ちている）
> dyn.load("libgcc_s_dw2-1.dll ")
として R 上に読み込む必要があるかもしれない．

Fortran で書かれたサブルーチンは，write 文など外しておく必要があるものがいくつ
かある．
適用例 (3 重対角行列の逆)
予め，tridiaginv.for を用意し，上記の要領で，tridiaginv.dll を作成しておく．
この二つのファイルは
http://www.st.hirosaki-u.ac.jp/~sugimoto/R_sample/tridiaginv.for
http://www.st.hirosaki-u.ac.jp/~sugimoto/R_sample/tridiaginv.dll
にアップしておきました．
> np=5 #3 重対角行列の次元
> A=1:5; B=cos(1:(np-1)/10) # 3 重対角行列の成分: A=対角成分 B=(i,i-1)成分ベクトル
> V=matrix(0,ncol=np,nrow=np) # 逆行列
> dyn.load("tridiaginv.dll")
> .Fortran("tridiaginv",as.double(A),as.double(B),as.integer(np),as.double(V))
[[1]]
[1] 1 2 3 4 5
[[2]]
[1] 0.9950042 0.9800666 0.9553365 0.9210610
[[3]]
[1] 5
[[4]]
[1] 2.49514383 -1.50265083 0.53325270 -0.13300051 0.02450032 -1.50265083
1.51019551 -0.53593012 0.13366830 -0.02462333 0.53325270
[12] -0.53593012 0.55228042 -0.13774628 0.02537455 -0.13300051 0.13366830
-0.13774628 0.29542996 -0.05442180 0.02450032 -0.02462333
[23] 0.02537455 -0.05442180 0.21002516
> V=.Fortran("tridiaginv",as.double(A),as.double(B),as.integer(np),as.double(V))[[4]]
上の V が逆行列の結果になる．計算時間は R と比べ，格段に速い．
2013.2.16
杉本知之