BugTrack-Ｒ備忘録/58

Ｒ備忘録 - 記事一覧

R（R言語）からCを呼ぶ（Linux編） - blas・lapack

• 投稿者： みゅ
• カテゴリ： なし
• 優先度： 普通
• 状態： 完了
• 日時： 2010年01月24日 20時53分37秒

内容

• blas・lapackを使ってみます．
• ソース → test.c(98)
• バグ報告や質問は別館でお願いします．

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

コメント