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

!!!LAPACKのQR分解を行う関数
*投稿者： みゅ
*カテゴリ： なし
*優先度： 普通
*状態： 完了
*日時： 2010年05月26日 15時03分32秒
//{{bugstate}}
!!内容
*自分用メモ
**ここで使用しているBLAS・LAPACK関数はATLAS化されています．
*ソース
 set.seed(1)
 nn <- 10000
 p <- 1000
 xx <- matrix(rnorm(nn*(p-1)), nn, p-1)
 dim(xx)
 xx1 <- cbind(1,xx)
 dim(xx1)
 b <- runif(p)
 yy <- cbind(1,xx) %*% b + rnorm(nn,0,0.2)
 dat <- data.frame(yy=yy, xx)
 dim(dat)

!!dgeqrf
*LWORKをN（列数）とおいた場合．
 system.time(
  res <- .Fortran("dgeqrf", M=as.integer(nn), N=as.integer(p),
              A=cbind(1,xx), LDA=as.integer(nn),
              TAU=double(p),
              WORK=double(p), LWORK=as.integer(p),
              INFO=integer(1))
  )
    ユーザ   システム       経過
     36.684      0.430     37.821
*LWORKを計算させる
 system.time({
  res <- .Fortran("dgeqrf", M=as.integer(nn), N=as.integer(p),
              A=cbind(1,xx), LDA=as.integer(nn),
              TAU=double(p),
              WORK=double(1), LWORK=-1L,  # LWORKを計算させる
              INFO=integer(1))
  print(res$WORK)
  res <- .Fortran("dgeqrf", M=as.integer(nn), N=as.integer(p),
              A=cbind(1,xx), LDA=as.integer(nn),
              TAU=double(p),
              WORK=double(res$WORK), LWORK=as.integer(res$WORK),
              INFO=integer(1))
  })
 [1] 36000
    ユーザ   システム       経過
     13.084      0.676      5.811

!!dgeqp3
*JPVT=1:Nと、おいた場合．
 system.time({
  res <- .Fortran("dgeqp3", M=as.integer(nn), N=as.integer(p),
              A=cbind(1,xx), LDA=as.integer(nn), JPVT=1L:p,
              TAU=double(p),
              WORK=double(1), LWORK=-1L,  # LWORKを計算させる
              INFO=integer(1))
  print(res$WORK)
  res <- .Fortran("dgeqp3", M=as.integer(nn), N=as.integer(p),
              A=cbind(1,xx), LDA=as.integer(nn), JPVT=1L:p,
              TAU=double(p),
              WORK=double(res$WORK), LWORK=as.integer(res$WORK),
              INFO=integer(1))
  })
 [1] 38036
    ユーザ   システム       経過
     13.040      0.556      5.621
*JPVT=０（integer(p)）と、おいた場合．
 system.time({
  res <- .Fortran("dgeqp3", M=as.integer(nn), N=as.integer(p),
              A=cbind(1,xx), LDA=as.integer(nn), JPVT=integer(p),
              TAU=double(p),
              WORK=double(1), LWORK=-1L,  # LWORKを計算させる
              INFO=integer(1))
  print(res$WORK)
  res <- .Fortran("dgeqp3", M=as.integer(nn), N=as.integer(p),
              A=cbind(1,xx), LDA=as.integer(nn), JPVT=integer(p),
              TAU=double(p),
              WORK=double(res$WORK), LWORK=as.integer(res$WORK),
              INFO=integer(1))
  })
 [1] 38036
    ユーザ   システム       経過
     18.213      0.576     14.622
**遅くなる
!!LINPACKのdqrdc2
 system.time(
  res <- .Fortran("dqrdc2", x=cbind(1,xx), ldx=as.integer(nn),
              n=as.integer(nn), p=as.integer(p),
              tol=1.0e-7, k=integer(1),
              qraux=double(p), jpvt=1L:p,
              work=double(p*2))
  )
    ユーザ   システム       経過
     35.556      0.282     35.833



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