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

!!!R（R言語）からCを呼ぶ（Linux編） - blas・lapack
*投稿者： みゅ
*カテゴリ： なし
*優先度： 普通
*状態： 完了
*日時： 2010年01月24日 20時53分37秒
//{{bugstate}}
!!内容
*blas・lapackを使ってみます．
*ソース → {{ref test.c}}
*バグ報告や質問は[別館|http://phoenixx.sakura.ne.jp/mt/R/2010/01/rc.html]でお願いします．
!!dgemv（blas）
*y := alpha*A*x + beta*yの計算を行う．
*計算した結果は「ｙ」に入ります．
*ここではalphaとbetaは「1.0」として実装しています．
*Cのソース
 //////////////////////////////////////////////////
 // ＢＬＡＳを使うテスト(dgemv)
 // y := alpha*A*x + beta*y
 long int gl_incx = 1;
 long int gl_incy = 1;
 double alphaOne = 1.0;
 double betaOne = 1.0;
 double betaZero = 0.0;
 char TRANSN = 'N';
 char TRANST = 'T';
 
 void _dgemv(long int m, long int n, double *A, long int lda, double *x, double *ans){
   dgemv_(&TRANSN, &m, &n, &alphaOne, A, &lda, x, &gl_incx, &betaOne, ans, &gl_incy);
 }
 
 SEXP Rdgemv01(SEXP ex_A, SEXP ex_x, SEXP ex_y)
 {
   SEXP ans;
   int m = INTEGER(GET_DIM(ex_A))[0];
   int n = INTEGER(GET_DIM(ex_A))[1];
   double *A;
   double *x;
   double *y;
   int lda = m;
   long int i;
   
   A = REAL(ex_A);
   x = REAL(ex_x);
   y = REAL(ex_y);
   PROTECT(ans = allocVector(REALSXP, m));
   for( i=0; i<m; i++){
     REAL(ans)[i] = y[i];
   }
   //dgemv_(&TRANSN, &m, &n, &alphaOne, A, &lda, x, &gl_incx, &betaOne, REAL(ans), &gl_incy);
   _dgemv(m, n, A, lda, x, REAL(ans));
   UNPROTECT(1);
   return(ans);
 }
*Rで実行
 > dyn.load("test.so")
 > is.loaded("Rdgemv01")
 [1] TRUE
 > set.seed(0)
 > A <- matrix(rnorm(12), 3, 4)
 > x <- rnorm(4)
 > y <- rep(0,3)
 > A %*% x + y
             [,1]
 [1,] -2.52945736
 [2,]  0.02833907
 [3,] -0.74987050
 > .Call("Rdgemv01", A, x, y)
 [1] -2.52945736  0.02833907 -0.74987050  //同じ値
!転置も考慮
*Aを転置して計算したい場合、Rで転置させることも可能ですが、Rで転置させるととても遅いので、blasにやってもらったほうがいいです．
 //////////////////////////////////////////////////
 // ＢＬＡＳを使うテスト(dgemv)
 // y := alpha*A*x + beta*y
 long int gl_incx = 1;
 long int gl_incy = 1;
 double alphaOne = 1.0;
 double betaOne = 1.0;
 double betaZero = 0.0;
 char TRANSN = 'N';
 char TRANST = 'T';
 
 //転置させないで計算
 void _dgemv(long int m, long int n, double *A, long int lda, double *x, double *ans){
   dgemv_(&TRANSN, &m, &n, &alphaOne, A, &lda, x, &gl_incx, &betaOne, ans, &gl_incy);
 }
 //転置して計算
 void _dgemvT(long int m, long int n, double *A, long int lda, double *x, double *ans){
   dgemv_(&TRANST, &m, &n, &alphaOne, A, &lda, x, &gl_incx, &betaOne, ans, &gl_incy);
 }
 
 SEXP Rdgemv01(SEXP ex_A, SEXP ex_x, SEXP ex_y, SEXP TRANS, SEXP NN)
 {
   SEXP ans;
   int m = INTEGER(GET_DIM(ex_A))[0];
   int n = INTEGER(GET_DIM(ex_A))[1];
   int nn = INTEGER(NN)[0];
   double *A;
   double *x;
   double *y;
   int lda = m;
   long int i;
   int trans = 0;
   
   if( LOGICAL(TRANS)[0] ){
     trans = 1;
   }
   A = REAL(ex_A);
   x = REAL(ex_x);
   y = REAL(ex_y);
   int mn = (trans==0)?m:n;
   PROTECT(ans = allocVector(REALSXP, mn));
   for( i=0; i<nn; i++){
     dcopy_(&mn, y, &gl_incx, REAL(ans), &gl_incy);
     if( trans == 0 ){
       _dgemv(m, n, A, lda, x, REAL(ans));
     } else{
       _dgemvT(m, n, A, lda, x, REAL(ans));
     }
   }
   UNPROTECT(1);
   return(ans);
 }
*Rの結果
 > dyn.load("test.so")
 > is.loaded("Rdgemv01")
 [1] TRUE
 >
 > set.seed(0)
 > A <- matrix(rnorm(12), 3, 4)
 > x <- rnorm(4)
 > y <- rep(1,3)
 > A %*% x + y
            [,1]
 [1,] -1.5294574
 [2,]  1.0283391
 [3,]  0.2501295
 > .Call("Rdgemv01", A, x, y, F, as.integer(1))
 [1] -1.5294574  1.0283391  0.2501295
 >
 > x1 <- rnorm(3)
 > y1 <- rep(1,4)
 > t(A) %*% x1 + y1
           [,1]
 [1,] 2.1888924
 [2,] 0.2801785
 [3,] 1.0261484
 [4,] 0.5773298
 > .Call("Rdgemv01", A, x1, y1, T, as.integer(1))
 [1] 2.1888924 0.2801785 1.0261484 0.5773298
 >
 >
 > system.time(for(i in 1:100000) A %*% x + y )
    ユーザ   システム       経過
       0.58       0.00       0.58
 > system.time(.Call("Rdgemv01", A, x, y, F, as.integer(100000)) )
    ユーザ   システム       経過
       0.06       0.00       0.06
 > system.time(for(i in 1:100000) t(A) %*% x1 + y1 )
    ユーザ   システム       経過
       2.51       0.00       2.53                  # 転置はとてもコストがかかる
 > system.time(.Call("Rdgemv01", A, x1, y1, T, as.integer(100000)) )
    ユーザ   システム       経過
       0.07       0.00       0.07
!!行列式を求めたい
*Cのソース
 //////////////////////////////////////////////////
 // 行列式のテスト
 // コレスキー分解を使用
 SEXP detTest01(SEXP Ra){
   char *UPLO = "U";
   int i;
   int n = INTEGER(GET_DIM(Ra))[0];
   int lda = n;
   SEXP ans;
   SEXP Rw;
   SEXP Rinfo;
   SEXP Rdet;
   PROTECT(Rw = allocMatrix(REALSXP, n, n));
   PROTECT(Rinfo = allocVector(INTSXP, 1));
   PROTECT(Rdet = allocVector(REALSXP, 1));
   PROTECT(ans = allocVector(VECSXP, 3));
   SET_VECTOR_ELT(ans, 0, Rw);
   SET_VECTOR_ELT(ans, 1, Rdet);
   SET_VECTOR_ELT(ans, 2, Rinfo);
   for( i=0; i<n*n; i++ ){
     REAL(Rw)[i] = REAL(Ra)[i];
   }
   dpotrf_(UPLO, &n, REAL(Rw), &lda, INTEGER(Rinfo));
   REAL(Rdet)[0] = 1.0;
   for( i=0; i<n; i++ ){
     REAL(Rdet)[0] = REAL(Rdet)[0] * REAL(Rw)[i+i*n]*REAL(Rw)[i+i*n];
   }
   UNPROTECT(4);
   return ans;
 }
*Rで実行
 > dyn.load("test.so")
 > A <- cov(matrix(rnorm(1000),200,5))
 > det(A)
 [1] 0.8832757
 > ret <- .Call("detTest01", A)
 > ret
 [[1]]
             [,1]        [,2]        [,3]        [,4]        [,5]
 [1,]  0.91836184  0.02698488  0.11704126  0.17378581 -0.05802576
 [2,]  0.02478189  1.01322217  0.01466400  0.08903543 -0.08439662
 [3,]  0.10748622  0.01801624  1.03501661 -0.08813723 -0.07886582
 [4,]  0.15959826  0.09490226 -0.06957777  0.91149304 -0.11018057
 [5,] -0.05328865 -0.08707835 -0.08965644 -0.11107615  1.07060425
 
 [[2]]
 [1] 0.8832757         # これが行列式
 
 [[3]]
 [1] 0



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