//このページが表示される方は、URLから「action=SOURCE&」を削除してみてください
[[Ｒ備忘録 - 記事一覧]]

!!!R（R言語）でロードできるライブラリ（.soファイル）のコンパイルのお作法？（Linux）
*投稿者： みゅ
*カテゴリ： なし
*優先度： 普通
*状態： 完了
*日時： 2009年06月03日 08時57分46秒
//{{bugstate}}
!!内容
*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
*R
 > .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

!!コメント
//{{comment}}