トップ 差分 一覧 ソース 検索 ヘルプ PDF RSS ログイン

RからBLAS/LAPACKを使う

RからC言語を使い、R内部の関数のみで作るよりも高速なプログラムを作成する方法を解説しています。
RjpWiki: Rから他言語利用 の応用編です。

最終更新時間:2007年05月31日 15時06分07秒


なぜRからCをつかうのか

Rの汎用性とCのスピードを同時に実現するためです。数値計算だけなら、MATLABでやるのが文法も分かりやすく手軽ですが、スピードと汎用性に欠けます。Rだけでやると、スピードはMATLABよりもだいぶ劣ります。Oxというのもあり、スピードが速いそうですが、汎用性に欠けます。特に、無償版はインタープリタ型ではないのと、グラフィックスと結果のデータの扱いに難があります。

別に、Fortranでもよいのですが、少々古く、最近はもっぱら数値計算用のみにつかわれているので、一般性を考えれば、Cのほうが無難だと思います。

BLASとLAPACKについて

Cでベクトル・行列演算を行う場合、BLAS/LAPACKという、Fortranで書かれたライブラリを使います。BLAS/LAPACKはRに実装されているので、RからC言語を使う場合、新たにBLAS/LAPACKをインストールしなくてもつかうことができます。

BLASとは、Fortranで書かれたベクトル・行列演算ライブラリです。LAPACKもFortranで書かれたライブラリで、最小二乗法や固有値問題ルーチンなどが実装されています。

C:\Program Files\R\rw2010\include\R_ext\BLAS.h


C:\Program Files\R\rw2010\include\R_ext\LAPACK.h

あたりを参考にしてください。BLAS/LAPACK関連のものは、すべてNetlibにあります。

(もしかしたら、R自身がBLAS/LAPACKを内装しているのだから、Rの関数をCの内部に組み込んでコンパイルするのとスピードはほとんど変わらない可能性はあります。ただ、実際のところは比較して見ないと分からないのと、BLAS/LAPACKを直接つかったほうが少なくとも遅くはならないというのは確かです。)

最初にすべきこと

Windows の場合は、以下の工程を踏む必要があります。

次に、BLASとLAPACKを使うため、"Makevars"という名前のファイルを作り、

PKG_LIBS = $(LAPACK_LIBS) $(BLAS_LIBS) $(FLIBS)

と記述して、Rのカレント・ディテクトリにおきます。

例題

例題1: 内積(1)

1) "myfunc.c"という名前のCファイルを作成します。

#include <R_ext/blas.h>
void myfunc(double *a, double *b, int *c, double *d)
{
 	int one=1;
	*d = F77_CALL(ddot)(c,a,&one,b,&one);
}

2) コマンドプロンプトから、

Rcmd SHLIB myfunc.c

と打ち込むことで、コンパイルできます。

3) Rから、

dyn.load("myfunc.dll")
myfunc <- function(a,b){
    .C("myfunc",
      arg1=as.double(a),
      arg2=as.double(b),
      arg3=as.integer( (length(a)) ), 
      arg4=as.double(0)
      )$arg4
 }と打ち込み、

例えば、

vec1<-c(1,2,3,4,5,6)
vec2<-c(2,3,4,5,6,8)
myfunc(vec1,vec2)

と打ち込めば、各々結果が返されます。

myfunc(c(1,2,3,4,5,6), c(2,3,4,5,6,8))

でも大丈夫です。

ちなみに、共用ライブラリを切り離すには、

dyn.unload("myfunc.dll")

と打ち込めばよい。

例題2: 内積(2)

"SEXP"を使う場合。

1) "myfunc2.c"という名前のCファイルを作成します。

#include <R.h>
#include <Rinternals.h>
#include <R_ext/blas.h>
SEXP myfunc2(SEXP a, SEXP b)
{
	int na, nb;
	na = length(a); nb = length(b);
	const int one = 1;
	SEXP ans;
	PROTECT(ans = allocVector(REALSXP, one));
	REAL(ans)[0] = F77_CALL(ddot)(&na,REAL(a),&one,REAL(b),&one);
	UNPROTECT(1);
	return(ans);
}

2) コマンドプロンプトから、

Rcmd SHLIB myfunc2.c

と打ち込むことで、コンパイルできます。

3) Rから、

dyn.load("myfunc2.dll")
myfunc2 <- function(a, b) .Call("myfunc2", a, b)

と打ち込み、例えば、

myfunc2(c(1,2,3),c(4,5,6))

とすれば結果が返されます。

C++の場合

C++で関数を書いた場合には、

extern "C"{}

でくくっておけばRから呼び出せます。

//myfunc2.cpp(C++独自の機能を使っていないのでよい例ではありません。)

extern "C"{
#include <R.h>
#include <Rinternals.h>
#include <R_ext/blas.h>
}

extern "C"{
SEXP myfunc2(SEXP a, SEXP b)
{
	int na, nb;
	na = length(a); nb = length(b);
	const int one = 1;
	SEXP ans;
	PROTECT(ans = allocVector(REALSXP, one));
	REAL(ans)[0] = F77_CALL(ddot)(&na,REAL(a),&one,REAL(b),&one);
	UNPROTECT(1);
	return(ans);
}
}

optim()について

optim()はRに組み込まれている汎用最適化関数で、多変数関数の最小値(マイナスをつければ最大値)が求められます。

RjpWikiの関数の最大・最小化の項を参照してください。BLAS/LAPACKの場合と同じく、Cで書かれたのプログラムの中に直接optim()のコード(正確には、そこから呼び出されるvmmin()、nmmin()やsamin())を組み込むことで、大規模なプログラムを組む場合の高速化が実現できます。

Nelder-Mead法のnmmin()をつかう場合

まず、"myfunc3.c"というファイル名で以下を保存します。

#include <R.h>
#include <Rinternals.h>
#include <R_ext/Applic.h>

double f(int n, double *x, void *dummy){
	double ans = 0;
	ans = x[0]*x[0] + (x[1]+1)*(x[1]+1) + 2;
	return ans;
}

SEXP myfunc3(SEXP dum)
{
	double x[2]={0,0};
	double opar[2]={0,0};
	double Fmin;
	int fncount, ifail;
  
	nmmin(2, x, opar, &Fmin, f,&ifail, 1.0e-8, 1.0e-8, NULL,
	1, 0.5, 2, 0,&fncount, 1000);
	
	SEXP ans;
	PROTECT(ans = allocVector(REALSXP, 3));
	REAL(ans)[0] = Fmin;
	REAL(ans)[1] = x[0];
	REAL(ans)[2] = x[1];
	UNPROTECT(1);
	return(ans);
}

次に、コマンドプロンプトから

Rcmd shlib myfunc3.c

でコンパイル。

結果はRのコンソールから

dyn.load("myfunc3.dll")
res3 <- .Call("myfunc3")
res3

と打ち込めば出ます。