数値解析の速度は申し分がないFortranですが、プロットなどのI/O周りの使い勝手は悪いです。Fortranでサブルーチンを書いて、Rから呼び出すのが良いかなとは思うのですが、Fortranの計算結果をファイルに書き出して、Rに読ますのもひとつの方法です。テキストで書き出してもよいのですが、効率の問題でバイナリーの方が良いでしょう。
Fortranからの出力
試しに以下のFortranの出力を、Rに読ませます。3×3の行列です。
program makeBin implicit none double precision, dimension(3, 3) :: m integer i, j do i = 1, 3 do j = 1, 3 m(i, j) = i + 3*(j - 1) end do end do open(17,file='mydata_seq.dat', form='unformatted', status='replace') write(17) m close(17) open(17,file='mydata_str.dat', form='unformatted', access = "stream", status='replace') write(17) m close(17) end program makeBin
2つファイルが作成されます。mydata_seq.datが伝統的なバイナリ・ファイルの出力のファイルで、mydata_str.datがFortran 2003以降で使えるC言語と同様のフォーマットのファイルです。
伝統的な方は access = "sequential"が省略されています。このsequential方式だと、writeごとにrecord markerをデータの前後に4バイトづつ付けて来ます。何かの役に立っているのかは分かりません。
C言語互換/stream形式の読み方
変数の型と行列の縦横サイズはあらかじめ把握しておく必要があります。ここではFortranでdouble precisionであったので、numeric型として読み込んでいます*1。
size <- c(3, 3) fname <- "mydata_str.dat" fs <- file.info(fname)[, "size"] is <- file(fname, "rb", blocking = FALSE) buf <- readBin(is, "raw", fs) v <- readBin(buf, numeric(), n = prod(size)); close(is) m <- matrix(v, size[1], size[2])
Fortranの方でaccess = "stream"が使えるのであれば、これが無難です。
mmapを使うこともできます。こちらの方が高速なので、読み込むデータが多い場合は重宝するかも知れません。
library(mmap) s <- mmap(fname, mode = double()) m <- matrix(s[], size[1], size[2]) munmap(s)
従来方式/sequential形式の読み方
無闇にwriteされていたら収拾がつかない気がしますが、1回であればファイルの冒頭と末尾を落とせば済みます。
size <- c(3, 3) fname <- "mydata_seq.dat" fs <- file.info(fname)[, "size"] is <- file(fname, "rb", blocking = FALSE) buf <- readBin(is, "raw", fs) buf2 <- buf[-c(1:4, (fs - 3):fs)] v <- readBin(buf2, numeric(), n = prod(size)); close(is) m <- matrix(v, size[1], size[2])
mmapもオプションを指定すれば使うことができます。
library(mmap) fname <- "mydata_seq.dat" fs <- file.info(fname)[, "size"] s <- mmap(fname, mode = double(), off = 4, len = fs - 8) m <- matrix(s[], size[1], size[2]) munmap(s)
変数を並べて出力した場合
Fortranでwrite(7) v1, v2というように変数を並べて出力した場合、v1の要素が出力されたあとに続けてv2の要素が出力されます。以下が単純な例になりますが、
program mmap implicit none integer, dimension(100) :: v1 integer i double precision, dimension(100) :: v2 do i = 1, 100 v2(i) = i v1(i) = 0.123 end do open(17,file='mmap_str.dat', form='unformatted', access = "stream", status='replace') write(17) v1, v2 close(17) end program mmap
R側で開始位置（offset）と長さ（len）を調整して、2つのベクトルとして扱うことになります。
fname <- "mmap_str.dat" library(mmap) v1 <- mmap(fname, mode = integer(), len = 400) # integerのサイズは4バイト v2 <- mmap(fname, mode = double(), off = 400, len = 800) # doubleのサイズは8バイト # munmap(s)
派生型（構造体）で出力した場合
以下のような派生型を使った場合は、Cの構造体を使った場合と同様にmmapで読み込めます。
program makeBinS implicit none integer i type :: sd integer :: a double precision :: b end type sd type(sd), dimension(100) :: s do i = 1, 100 s(i)%a = i s(i)%b = 1.23 end do open(17,file='mydata_sd.dat', form='unformatted', access = "stream", status='replace') write(17) s close(17) end program makeBinS
mmapを使わないと、int型の読み込みreadBinとdouble型読み込みのreadBinを2つ並べてループするようになって煩雑になると思います。
fname <- "mydata_sd.dat" library(mmap) record.type <- struct(a = integer(), b = double()) s <- mmap(fname, mode = record.type) # munmap(s)
