プロセッサのマルチコア化が進んだ現代なので、時間のかかる計算では並列処理をするコードを書くのが望ましいです。スレッドやセマフォの制御を直接プログラミングすると骨が折れるのですが、シングルコア向け逐次処理コードを並列処理に変換してくれるOpenMPという仕組みがあります。Cで書いたRの拡張で使ってみましょう。
1. Makevarsの編集
gccでOpenMPを使うためには、引数に-fopenmpフラグを加える必要があります。R CMDでこの引数をつけるためには、LinuxやMacOS Xでは$HOME/.R/Makevarsに、MS-Windowsでは%HOME%*1/.R/Makevars.winと言うテキストファイルに、
CFLAGS= -fopenmp
CXXFLAGS= -fopenmp
# PKG_CXXFLAGSか、PKG_CXXFLAGSに-std=c++11と書いてCXX11FLAGSに-fopenmpをつける（C++14/17でも同様）ほうがよいときもあります。
と言うようにかいておきます。なお、-O3オプションなども同様につきます。
2. Cのソースコードの編集
ループされていて一定以上の処理時間がかかるコードの例として、2つの一行を除けば同じ（無意味な）関数をOMPexample.cとして保存します。
#include<R.h>
#include<Rinternals.h>SEXP example(SEXP m){
int i, j, n;
double s1 = 0, s2 = 0;
int nor, noc;
SEXP rv, dim;/* 第一引数が行列か確認 */
if(!isMatrix(m))
error("A matrix is required for the first argument.");/* 行列の行数、列数を得る */
dim = getAttrib(m, R_DimSymbol);
nor = INTEGER(dim)[0]; /* 行数 */
noc = INTEGER(dim)[1]; /* 列数 */for(j=0;j<noc;j++){
for(i=0;i<nor;i++){
s1 += REAL(m)[i+j*nor]*(i+1)/(j+1); /* 行列と言っても一次元配列 */
s2 -= REAL(m)[i+j*nor]*(i+1)/(j+1);
}
}rv = PROTECT(allocVector(REALSXP, 2));
REAL(rv)[0] = s1;
REAL(rv)[1] = s2;
UNPROTECT(1);return rv;
}SEXP OMPexample(SEXP m){
int i, j, n;
double s1 = 0, s2 = 0;
int nor, noc;
SEXP rv, dim;/* 第一引数が行列か確認 */
if(!isMatrix(m))
error("A matrix is required for the first argument.");/* 行列の行数、列数を得る */
dim = getAttrib(m, R_DimSymbol);
nor = INTEGER(dim)[0]; /* 行数 */
noc = INTEGER(dim)[1]; /* 列数 */#pragma omp parallel for private(i, j) reduction(+:s1) reduction(-:s2)
for(j=0;j<noc;j++){
for(i=0;i<nor;i++){
s1 += REAL(m)[i+j*nor]*(i+1)/(j+1); /* 行列と言っても一次元配列 */
s2 -= REAL(m)[i+j*nor]*(i+1)/(j+1);
}
}rv = PROTECT(allocVector(REALSXP, 2));
REAL(rv)[0] = s1;
REAL(rv)[1] = s2;
UNPROTECT(1);return rv;
}
大学の情報基盤センターなるところが書いたネット上にあるOpenMPの解説を読めばすぐ分かりますが、プリプロセッサ命令
#pragma omp parallel for private(i, j) reduction(+:s1) reduction(-:s2)
がOpenMPの記述になります。private(i, j)は並列処理で変数iとjは、それぞれ異なる値をシェアしない変数として扱うことを意味し、reduction(+:s1)は、それぞれの並列処理で加算されていくそれぞれの変数s1は、並列処理の終了時に一つに合算されることを意味します。reduction(-:s2)は、加算ではなく減算処理をあわせます。
なお、OpenMPが有効でない場合は、このプリプロセッサは無視されて逐次処理されるコードが出力されます。
3. Cのソースコードをコンパイルする
OMPexample.cを用意したら、以下のようにすると、OMPexample.so（もしくはOMPexample.dll）と言うRの拡張が出来上がります。
R CMD SHLIB OMPexample.c
Rtoolsが入っていないと動かないので注意してください。
4. Rから拡張を呼び出して使う
OpenMPを使ったOMPexample関数と、使っていないexample関数の速度比較をRでしてみましょう。
dll <- dyn.load(paste("OMPexample", .Platform$dynlib.ext, sep = ""))
example <- function(m) .Call("example", m)
OMPexample <- function(m) .Call("OMPexample", m)set.seed(1002)
n <- 5000 # 一瞬で処理が終わる人は増やすと時間をかけられます
m <- matrix(rnorm(n*n), n, n)print(system.time({
cat("No OpenMP: ", example(m), "\n")
}))print(system.time({
cat("OpenMP: ", OMPexample(m), "\n")
}))
搭載されているプロセッサ数に応じて速くなります。並列化と言っても、手軽に試せますね。
A. Rcpp/Eigen
Rcppの行列はスレッドセーフではないのでOpenMPの並行処理の中で使えませんが、Eigenは最初にEigen::initParallel();をして行列をスレッドセーフにしておくだけで、OpenMPの並列処理の中で行列を使うことができます。つまり、Rcppの行列からEigenの行列を作成して処理を行うことで、RcppからOpenMPを使うことができます。RcppPararellを使うべきと思うかも知れませんが、RcppPararellはWorkerクラスを定義しないといけないので記述が煩雑な一方、Eigen/OpenMPと比較して速度的なアドバンテージはほぼ無いようです。
