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

BugTrack-Ｒ備忘録/28

R（R言語）でロードできるライブラリ（.soファイル）のコンパイルのお作法？（Linux）

投稿者：みゅ
カテゴリ：なし
優先度：普通
状態：完了
日時： 2009年06月03日 08時57分46秒

内容

R_HOMEにあるMakefileをインクルードして、Makeすると、便利

Makefile雛形

R_SHARE_DIR = ${R_HOME}/share/
R_INCLUDE_DIR = ${R_HOME}/include

#include $(R_HOME)/etc${R_ARCH}/Makeconf
include ${R_HOME}/etc${R_ARCH}/Makeconf

TARGET = kawase007.so

COBJS = kawase007_main.o kawase007.o Rconnection.o
#HEADER = 
#BASIS = *.h

$(TARGET): $(COBJS)
	g++ -std=gnu99 $(SHLIB_LDFLAGS) $(LDFLAGS) -o $@ $(COBJS) $(ALL_LIBS)  -llapack -lg2c -lm
#	$(SHLIB_LINK) -o $@ $(COBJS) $(ALL_LIBS)  -llapack -lg2c -lm
# 通常なら上記コメントアウトしている行で、コンパイルできるはず
# $(SHLIB_LINK)で指定されているリンカー（gcc）を使うとエラーが出る場合は、変更してみるとよい
# c++などが、混ざっている場合など？
# 上記ではg++に変更している

clean :
	rm -f *.o

blasとか、lapack関連のライブラリを作成したい場合

blaslapack.so : blaslapack.o
	$(SHLIB_LINK) -o $@ blaslapack.o $(ALL_LIBS)  -llapack -lg2c -lm

atlasとかだと

blaslapack.so : blaslapack.o
        $(SHLIB_LINK) -o $@ blaslapack.o $(ALL_LIBS) -L/usr/local/atlas/lib -llapack -lptcblas -lptf77blas -latlas

${R_HOME}/etc${R_ARCH}/Makeconfを、もっと詳しくみると、こんなことをわざわざ、書く必要もないかもしれないが．

2009/10/24追加

Gentooでの例①

test01.c

#include <stdio.h>
#include "R_ext/BLAS.h"
void test01(double a[]){
        a[0] = a[0] * 2;
        return;
}

void test02(int *n, double *a, double x[], double y[]){
        int incx = 1;
        int incy = 1;
        daxpy_(n, a, x, &incx, y, &incy);
}

Makefile

R_SHARE_DIR = ${R_HOME}/share/
R_INCLUDE_DIR = ${R_HOME}/include

#include $(R_HOME)/etc${R_ARCH}/Makeconf
include ${R_HOME}/etc${R_ARCH}/Makeconf

%.so : %.o
        $(SHLIB_LINK) -o $@ $< $(ALL_LIBS)

clean :
        rm -f *.o

make

$ make test01.so
gcc -std=gnu99 -I/usr/local/lib/R/include -I/usr/local/include -fpic -g -O2 -c test01.c -o test01.o
gcc -std=gnu99 -shared -L/usr/local/lib -o test01.so test01.o   -L/usr/local/lib/R/lib -lR

> .C("test01", 3)
[[1]]
[1] 6
> .C("test01", a=c(3,9))
$a
[1] 6 9
> .C("test02", n=as.integer(3), a=0.1, x=c(3,6,9), y=c(0,0,0))
$a
[1] 6 9
> .C("test02", n=as.integer(3), a=0.1, x=c(3,6,9), y=c(0,0,0))
$n
[1] 3

$a
[1] 0.1

$x
[1] 3 6 9

$y
[1] 0.3 0.6 0.9

Gentooでの例②

test02.c

#include <stdio.h>
#include "R_ext/BLAS.h"
void test01(double a[]){
        a[0] = a[0] * 2;
        return;
}

void test02(int *n, double *a, double x[], double y[]){
        int incx = 1;
        int incy = 1;
        daxpy_(n, a, x, &incx, y, &incy);
}

void test03(int *n, double a[], int *info){
        char *UPLO = "U";
        int lda = *n;
        dpotrf_(UPLO, n, a, &lda, info);
}

void test04(int *n, double a[], int *info){
        char *JOBZ = "N";
        char *UPLO = "U";
        int lda = *n;
        double w[*n];
        int lwork = 3*(*n)-1;
        double work[lwork];
        dsyev_(JOBZ, UPLO, n, a, &lda, w, work, &lwork, info);
}

実行結果

> .C("test04", n=as.integer(dim(A)[1]), a=A, w=double(dim(A)[1]), info=integer(1))
$n
[1] 5

$a
             [,1]        [,2]         [,3]        [,4]       [,5]
[1,]  0.914646112  0.11281429 -0.122770717  0.07373369 -0.0662105
[2,] -0.088367911  0.89860486 -0.135535733  0.39746774  0.6268844
[3,]  0.084026524 -0.03413390  0.874546126 -0.19257233  0.3982850
[4,] -0.012745596  0.21714071  0.008365559  0.94047705  0.1168292
[5,]  0.009942577 -0.09413683 -0.059808934 -0.03333698  1.1367381

$w
[1] 0.6459856 0.8160620 0.9824064 1.0956348 1.2249235

$info
[1] 0

> eigen(A)
$values
[1] 1.2249235 1.0956348 0.9824064 0.8160620 0.6459856

$vectors
            [,1]       [,2]      [,3]        [,4]       [,5]
[1,]  0.22305982  0.5677500 0.6318367 -0.45577070  0.1447748
[2,] -0.48195385 -0.3237559 0.2476938 -0.04993972  0.7739874
[3,] -0.03343607  0.3956827 0.2586180  0.87261401  0.1182319
[4,] -0.43461646 -0.3182952 0.5976240  0.03090832 -0.5930320
[5,]  0.72660537 -0.5612178 0.3396947  0.16543443  0.1196585

なぜ、リンク時にblasとかlapackとかリンクしてないのに、使えるのかと思ったら・・・

# ldd /usr/local/lib/R/lib/libR.so
       linux-gate.so.1 =>  (0xffffe000)
       libblas.so.0 => /usr/lib/libblas.so.0 (0xb7c03000)
       liblapack.so.0 => /usr/lib/liblapack.so.0 (0xb7780000)
       libatlas.so.0 => /usr/lib/libatlas.so.0 (0xb73f0000)
       libgfortran.so.1 => /usr/lib/gcc/i486-pc-linux-gnu/4.1.2/libgfortran.so.1 (0xb7375000)
       libm.so.6 => /lib/libm.so.6 (0xb734e000)
       libreadline.so.5 => /lib/libreadline.so.5 (0xb731b000)
       libdl.so.2 => /lib/libdl.so.2 (0xb7317000)
       libc.so.6 => /lib/libc.so.6 (0xb71e7000)
       /lib/ld-linux.so.2 (0xb7f47000)
       libgcc_s.so.1 => /usr/lib/gcc/i486-pc-linux-gnu/4.3.2/libgcc_s.so.1 (0xb71d9000)
       libcblas.so.0 => /usr/lib/libcblas.so.0 (0xb71b9000)
       libncurses.so.5 => /lib/libncurses.so.5 (0xb7173000)

# nm /usr/local/lib/R/lib/libR.so | grep daxpy
         U daxpy_
# nm /usr/lib/liblapack.a | grep dpotrf
ATL_dpotrf.o:
00000000 T ATL_dpotrf
         U ATL_dpotrfL
         U ATL_dpotrfU
ATL_dpotrfL.o:
00000000 T ATL_dpotrfL
ATL_dpotrfU.o:
00000000 T ATL_dpotrfU
         U ATL_dpotrf
ATL_f77wrap_dpotrf.o:
         U ATL_dpotrf
00000000 T atl_f77wrap_dpotrf_
         U ATL_dpotrf
clapack_dpotrf.o:
         U ATL_dpotrf
00000000 T clapack_dpotrf
         U dpotrf_
dpotrf.o:
         U atl_f77wrap_dpotrf_
00000000 T dpotrf_
         U dpotrf_
         U dpotrf_
         U dpotrf_

ということは、今呼ばれている関数は、atlasでないってことか．

Gentooでの例③

test03.c

#include <stdio.h>
#include "Rdefines.h"
#include "R_ext/BLAS.h"
#include "R_ext/Lapack.h"

・・・

void test04(int *n, double a[], double w[], int *info){
        char *JOBZ = "N";
        char *UPLO = "U";
        int lda = *n;
        //double w[*n];
        int lwork = 3*(*n)-1;
        double work[lwork];
        dsyev_(JOBZ, UPLO, n, a, &lda, w, work, &lwork, info);
}

SEXP test05(SEXP Rn, SEXP Ra, SEXP Rw, SEXP Rinfo){
        char *JOBZ = "N";
        char *UPLO = "U";
        int n = INTEGER(Rn)[0];
        int lda = n;
        int lwork = 3*(n)-1;
        double work[lwork];
        dsyev_(JOBZ, UPLO, &n, REAL(Ra), &lda, REAL(Rw), work, &lwork, INTEGER(Rinfo));
        return Rw;
}

結果

> dyn.load("test03.so")
> .C("test04", n=as.integer(dim(A)[1]), a=A, w=double(dim(A)[1]), info=integer(1))
$n
[1] 5

$a
             [,1]        [,2]         [,3]        [,4]       [,5]
[1,]  0.914646112  0.11281429 -0.122770717  0.07373369 -0.0662105
[2,] -0.088367911  0.89860486 -0.135535733  0.39746774  0.6268844
[3,]  0.084026524 -0.03413390  0.874546126 -0.19257233  0.3982850
[4,] -0.012745596  0.21714071  0.008365559  0.94047705  0.1168292
[5,]  0.009942577 -0.09413683 -0.059808934 -0.03333698  1.1367381

$w
[1] 0.6459856 0.8160620 0.9824064 1.0956348 1.2249235

$info
[1] 0
> .Call("test05", n=as.integer(dim(A)[1]), a=A, w=double(dim(A)[1]), info=integer(1))
[1] 0.6459856 0.8160620 0.9824064 1.0956348 1.2249235
> eigen(A)
$values
[1] 1.0956954+0.1880367i 1.0956954-0.1880367i 0.9341475+0.1143243i
[4] 0.9341475-0.1143243i 0.7053265+0.0000000i

$vectors
                      [,1]                  [,2]                    [,3]
[1,] -0.5096168+0.1154550i -0.5096168-0.1154550i  0.79037907+0.00000000i
[2,] -0.1020337-0.4360399i -0.1020337+0.4360399i  0.03818506+0.11768606i
[3,]  0.5414950+0.0000000i  0.5414950+0.0000000i -0.06175406-0.58207539i
[4,] -0.2286156-0.3090191i -0.2286156+0.3090191i  0.07552457-0.01562969i
[5,]  0.2889004+0.0445094i  0.2889004-0.0445094i  0.03088052-0.10230174i
                        [,4]           [,5]
[1,]  0.79037907+0.00000000i  0.36196003+0i
[2,]  0.03818506-0.11768606i -0.59582930+0i
[3,] -0.06175406+0.58207539i  0.43142099+0i
[4,]  0.07552457+0.01562969i  0.57154733+0i
[5,]  0.03088052+0.10230174i -0.03437971+0i

なんかめちゃくちゃあやしいんですが？
- ああ、しまった．こっちの「.Call」の呼び方だと、A自体を変更してしまうんだった．あぶないあぶない．

Gentooでの例④

#include <stdio.h>
#include "Rdefines.h"
#include "R_ext/BLAS.h"
#include "R_ext/Lapack.h"

・・・

SEXP test06(SEXP Ra){
  char *JOBZ = "N";
  char *UPLO = "U";
  int n = INTEGER(GET_DIM(Ra))[0];
  int lda = n;
  int lwork = 3*(n)-1;
  double work[lwork];
  SEXP Rw;
  int info;
  PROTECT(Rw = allocVector(REALSXP, n));
  dsyev_(JOBZ, UPLO, &n, REAL(Ra), &lda, REAL(Rw), work, &lwork, &info);
  UNPROTECT(1);
  return Rw;
}

実行結果

> dyn.load("test04.so")
> A <- cov(matrix(rnorm(500),100,5))
> .C("test04", n=as.integer(dim(A)[1]), a=A, w=double(dim(A)[1]), info=integer(1))
$n
[1] 5

$a
             [,1]       [,2]        [,3]       [,4]        [,5]
[1,]  0.834463115 0.15527427 -0.48993713 0.03185757  0.04401092
[2,] -0.008567303 0.90121558 -0.17016539 0.01605062  0.43986948
[3,]  0.136655892 0.09518474  0.97289541 0.22299300  0.22452943
[4,] -0.062643576 0.06514187  0.11326410 0.96736289 -0.18335848
[5,]  0.012954799 0.12947744  0.06609119 0.11099574  1.05222868

$w
[1] 0.6190685 0.7479007 0.9502250 1.1082119 1.3027596

$info
[1] 0

> .Call("test06", a=A)
[1] 0.6190685 0.7479007 0.9502250 1.1082119 1.3027596

Ｒ備忘録 /状態空間モデリング/donlp2/その他のメモ

BugTrack-Ｒ備忘録/28

R（R言語）でロードできるライブラリ（.soファイル）のコンパイルのお作法？（Linux）

内容

Makefile雛形

blasとか、lapack関連のライブラリを作成したい場合

atlasとかだと

2009/10/24追加

Gentooでの例①

Gentooでの例②

Gentooでの例③

Gentooでの例④

コメント

メニュー

サイト内検索

人気検索キーワード

人気の２０件

最新


		2024-4
日	月	火	水	木	金	土
	1	2	3	4	5	6
7	8	9	10	11	12	13
14	15	16	17	18	19	20
21	22	23	24	25	26	27
28	29	30