1954年に考案されたFortranは最古の高級言語かつ、数値解析では最速のプログラミング言語として知られています1が、全般的に使い勝手が悪い2です。Rは柔軟性がある3プログラミング言語で使い勝手に優れていますが、数値解析をさせるとFortranと比較して10倍から100倍は低速といった感じです4。数値解析と統計処理は浮動小数点演算を扱うと言う意味では似ているのですが、RとFortranのコンセプトは真逆に近いです5。これは連携がとれれば相乗効果を発揮することを意味します。そして、RからFortranのサブルーチンを呼び出すのは簡単です。
1 準備
MS-WindowsではRとRtoolsをインストールして、環境変数が適切に設定されているか確認してください。この文書の記述時点での最新のR
4.2.2とrtools
4.2の標準インストールであれば、環境変数RTOOLS42_HOMEがC:\rtools42になります。また、コマンドプロンプトからRを呼び出せるように、PATHに%R_HOME%\binを追加しておいてください6。
LinuxであればOSのパッケージマネージャーを用いるなどしてgfortranとRが動く環境にしてください。
インストールが済んだら、OpenMPが有効になるようにオプションを設定します。Windowsの場合は%R_USER%\.R\Makevars.winを作成します。Linuxの場合は~/.R/Makevarsです。Rcppなどで開発がしたことがある人にはお馴染みだと思います。
Windowsでこのファイルが存在しない場合は、Rで
directory <- file.path(Sys.getenv("R_USER"), ".R")
if(!dir.exists(directory)) dir.create(directory)
cat("CFLAGS = $(SHLIB_OPENMP_CFLAGS)\nFCFLAGS = $(SHLIB_OPENMP_FFLAGS)",
file=paste(directory, "Makevars.win", sep="/"))とすれば作られます。場所はprint(directory)をして確認してください。
作ったMakevars.winの中身のCFLAGSはCコンパイラに対するオプション、FCFLAGSはgFortran（のfree
source
formのコンパイル）に対するオプションになります。それぞれ$(SHLIB_OPENMP_CFLAGS)と$(SHLIB_OPENMP_FFLAGS)を指定しておくと、OpenMPが有効になります。Cコンパイラは不要に思えるかもですが、gFortranがつくったオブジェクト・ファイルをリンクするのに使います。
なお、Rの標準構成でも使っているBLASとLAPACKのサブルーチンは、FCFLAGSに特段の指定なくそのまま使えます。
2 プログラムの作成
設定が機能しているか確認していきましょう。
2.1 RとFortranの連携の確認
Fortranプログラマではないので、Calling Fortran subroutines from Rを少し改変して使います。
まず、以下の内容のfortsub.f90と言うファイルを作ります。
subroutine minus1(m,n,x,y)
integer,intent(in) :: m,n
double precision,intent(in) :: x(m,n)
double precision,intent(out) :: y(m,n)
y = x - 1
end subroutine minus1次に、以下の内容のfortsub.Rと言うファイルを作ります。
dll <- paste("fortsub", .Platform$dynlib.ext, sep="")
dyn.load(dll)
m <- 3
n <- 2
result <- .Fortran("minus1",
as.integer(m), # 型は明示しておかないと連携に失敗しがち
as.integer(n),
matrix(as.numeric(seq(1, m*n)),c(m,n)),
matrix(0, m, n))
dyn.unload(dll) # unloadしないとFortranのリコンパイル前にRを終了する必要が出てくる
print(result[[4]]) # 戻り値はリスト構造コマンドプロンプトで2つのファイルを置いたディレクトリに移動して、
R CMD SHLIB fortsub.f90
R --slave -f fortsub.Rとすれば実行でき、
[,1] [,2]
[1,] 0 3
[2,] 1 4
[3,] 2 5と結果を得ることができます。seq(1, m*n)はc(1,2,3,4,5,6)ですから、Fortranで-1できていることが分かります。連携成功です。
2.2 OpenMPが有効になっていることの確認
Fortranユーザーは速度を気にする数値解析ユーザーですから、やはり2022年時点ではOpenMPが必須でしょう。OpenMPが動いていることを確認します。OpenMP入門: 初めてのOpenMPプログラム - Hello World in OpenMPのソースコードを少し改変して試しましょう。
まず、以下の内容のhelloomp.f90と言うファイルを作ります。
subroutine helloomp
!$ use omp_lib
print *, "START"
!$omp parallel
print '("omp_get_num_threads = ", i2)', omp_get_num_threads()
print '("omp_get_thread_num = ", i2)', omp_get_thread_num()
!$omp end parallel
print *, "END"
end subroutine helloomp次に、次の内容のhelloomp.Rと言うファイルを作ります。
dll <- paste("helloomp.dll", .Platform$dynlib.ext, sep="")
fs <- dyn.load(dll)
result <- .Fortran("helloomp")
dyn.unload(dll)コマンドプロンプトで2つのファイルを置いたディレクトリに移動して、
R CMD SHLIB helloomp.f90
set OMP_NUM_THREADS=3
R --slave -f helloomp.Rとすれば実行でき、
START
omp_get_num_threads = 3
omp_get_thread_num = 0
omp_get_num_threads = 3
omp_get_thread_num = 2
omp_get_num_threads = 3
omp_get_thread_num = 1
ENDと結果を得ることができます。環境変数で3スレッドの利用を指定して、3種類のスレット番号の出力が出ているので、マルチスレッドになっていることが分かりますね。
なお、Rguiからhelloomp.Rを呼んでも、後述する理由で何も表示されません。
2.3 R経由でコンソール出力をする
Rguiから呼ばれた場合、Fortranのprintやwriteは画面に文字を表示しません。この問題を解決するために、R経由で表示するためのサブルーチンが用意されているので試してみましょう。
まず、以下の内容のファイルprint_via_R.f90を用意します。
subroutine printviar(x, n)
double precision,intent(in) :: x
integer,intent(in) :: n
call labelpr('Printing via R', -1) ! 文字列を表示
call dblepr1('x', 1, x) ! double型を表示
call intpr1("n", -1, n) ! integer型を表示
end subroutine printviarサブルーチンlabelprとdblepr1とintpr1の第1引数はラベルで、第2引数はラベルの長さです。-1を入れておくと文字列の長さを数えて全部表示してくれます。第3引数は中身を表示したい変数です。
次に、以下の内容のprint_via_R.Rを用意します。
dll <- paste("print_via_R", .Platform$dynlib.ext, sep="")
fs <- dyn.load(dll)
r <- .Fortran("printviar",
as.numeric(pi),
as.integer(1))
dyn.unload(dll)2つのファイルを置いたディレクトリで、print_via_R.f90をコンパイルした後にRguiを呼び出します。
R CMD SHLIB print_via_R.f90
R --slave -e "system(file.path(R.home(\"bin\"), \"Rgui.exe\"))"Rguiからprint_via_R.f90を呼ぶと、しっかり結果が戻ります。
source("print_via_R.R")Printing via R
x
[1] 3.141593
n
[1] 1なお、R
4.0以前はdobule型/real型/integer型の配列を表示するdblepr/realpr/intprだけが用意されており、R
4.0からラベルのみ表示するlabelprと、スカラーを表示できるdblepr1/realpr1/intpr1が用意されました。
2.4 Rでエラーと警告を出す
サブルーチンが提供されています。
call rwarn("This calcuration might be stopped!")
call rexit("An error occured!")2.5 ユーザー中断の許可
Fortranのコードの実行中はユーザー中断が許可されないのですが、ループ中にcall rchkusr()を挟んでおくと、CTRL+Cなどが押されていれば、そこで中断されます。
3 Fortran 90以降を利用するときの注意
?.Fortranをしたら書いてある注意事項ですが、小文字でアンダースコアなしの名前にした方が無難なようです。
Fortran
2003からiso_c_bindingでsubroutineにbind(c)をつけておくとCのオブジェクトファイルを生成するようになったので、“.C”関数でCの関数として呼び出すと言う手もあります。
SUBROUTINE do_exchange(a, b) BIND(C, NAME='exchange')
USE iso_c_binding
REAL(c_double), intent(inout) :: a, b
REAL(c_double) :: c
c = a
a = b
b = c
ENDと言う風にdo_exchange.f90をつくります。.f03は\R CMD SHLIBがフォートランのファイルとして認識しないので、.f90にしています。
dll <- "iso_c_binding.dll"
fs <- dyn.load(dll)
print(.C("exchange", as.double(pi), as.double(-1)))
dyn.unload(dll).CでBIND(C, NAME="...")したFortranのコードを呼ぶときは、NAMEで指定した名前で呼ぶことに注意してください。
4 RがCに提供しているAPIを呼び出す
Rオブジェクトを作るAPIの利用などは無理がありそうですが、算術関係の関数は便利なこともあるでしょう。2つ利用例をあげておきます。
4.1 引数の無いCの関数
Rのドキュメントに以下の乱数の生成のコードがあります。乱数生成は初期化の問題があるので、APIを使った方がよいでしょう。
module rngfuncs
use iso_c_binding
interface
double precision function unifRand() bind(C, name = "unif_rand")
end function unifRand
subroutine getRNGseed() bind(C, name = "GetRNGstate")
end subroutine getRNGseed
subroutine putRNGseed() bind(C, name = "PutRNGstate")
end subroutine putRNGseed
end interface
end module
subroutine make_a_number(x)
use rngfuncs
REAL(c_double) :: x
call getRNGseed()
y = unifRand()
call putRNGSeed()
end subroutineDLLの作り方とRからの呼び出し方は、これまでと同様です。ファイル名をmake_a_number.f90とでもして試してください。
4.2 引数のあるCの関数
非ポインター引数のあるCの関数を呼び出す場合、ややこしい記述になります。Fortranの関数呼び出しの引数は、参照渡し、Cで言うポインターになっています。ポインターを引数に持ち、目的の関数を内部で呼び出すCの中継関数を用意して、Fortranからは中継関数を呼び出すようにします。
以下の内容のrelay_pnorm.cを書きます。
#include<R.h>
#include<Rmath.h>
double relay_pnorm(double *x, double *mu, double *sigma, int *lt, int *lg){
return pnorm(*x, *mu, *sigma, *lt, *lg);
}以下の内容のcall_c_func.f90を書きます。
module interface_for_api
use iso_c_binding
interface
REAL(c_double) function pnorm(x, mu, sigma, lt, lg) bind(C, name = "relay_pnorm")
import c_double, c_int ! c_*とFortranの型の対応はiso_c_bindingを参照
implicit none
REAL(c_double), intent(in) :: x, mu, sigma
INTEGER(c_int), intent(in) :: lt, lg
end function
end interface
end module
subroutine ipnorm(d, x1, x2)
use interface_for_api
implicit none
REAL(c_double), intent(out) :: d
REAL(c_double), intent(in) :: x1, x2
REAL(c_double), parameter :: mu = 0.0, sigma = 1.0
INTEGER(c_int), parameter :: lt = -1, lg = 0 ! TRUEとFALSEに対応
d = pnorm(max(x1, x2), mu, sigma, lt, lg) - pnorm(min(x1, x2), mu, sigma, lt, lg)
end subroutine以下のようにDLLを作成します。
R CMD SHLIB call_api_for_c.f90 relay_pnorm.cRからの呼び出し方は、これまでと同様です。Rを起動して、以下を入力したら計算します。
dll <- paste("call_c_func", .Platform$dynlib.ext, sep="")
fs <- dyn.load(dll)
r <- .Fortran("ipnorm", numeric(1), as.numeric(-100), as.numeric(0))
print(r) # リストの最初が戻り値になる
dyn.unload(dll)5 まとめ
無制約とはいかないのですが、容易に連携できるので、Fortranで書いたサブルーチンをRに流用するのには使えると思います。
アルゴリズムを揃えた状態では（恐らく最適化のレベルが高いために）CやC++よりもパフォーマンスが高く、OpenMPを利用することで数値解析においてはマルチプロセッサにも容易に対応できます。↩︎
ファイル入出力やテキスト処理において最近人気のプログラミング言語に備わっているような正規表現のような機能がありませんし、そもそも設計が古いです。パッケージ配布システムも無く、埋め込みSQLが使えたりしますがデータベースとの接続も容易ではありません。プロットの機能もなく、数値解析ユーザーには不要でしょうが、対話型インターフェイスによる実行もできません。↩︎
簡潔だが十分に応用の利く言語仕様であること、容易にC/C++拡張が書けること、パッケージ配布システムが整っていることが使い勝手をよくしています。↩︎
線形代数演算ライブラリBLASに任せるベクトルや行列の処理、CやC++のコードに任せる拡張の処理は高速ですが、R自体で複雑なアルゴリズムを書くと関数型言語の限界により文字通り桁違いに遅くなります。↩︎
GUIのサポートが弱くアプリケーションの開発と配布では使われないと意味ではFortranとRは似ています。↩︎
%R_HOME%は標準インストール時はC:\Program Files\R\R-4.2.2\binになります。↩︎