BugTrack-Ｒ備忘録/65

Ｒ備忘録 - 記事一覧

LAPACKのQR分解を行う関数

• 投稿者： みゅ
• カテゴリ： なし
• 優先度： 普通
• 状態： 完了
• 日時： 2010年05月26日 15時03分32秒

内容

• 自分用メモ
  • ここで使用している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

コメント