MAGMA 和 Rcpp 用于 R 中的线性代数
Posted
技术标签:
【中文标题】MAGMA 和 Rcpp 用于 R 中的线性代数【英文标题】:MAGMA and Rcpp for linear algebra in R 【发布时间】:2013-08-23 11:43:30 【问题描述】:我想知道是否有人尝试使用Rcpp 和MAGMA 通过使用 CPU 和 GPU 来加速 R 中的线性代数运算?我上个月尝试了culatools,它与Rcpp (link) 一起使用,但是culatools 是一种商业产品,需要花钱才能访问所有功能。
【问题讨论】:
【参考方案1】:在修改了 culatools 之后,使用 Rcpp 和 MAGMA 非常简单。这是 .cpp 文件:
#include<Rcpp.h>
#include<magma.h>
using namespace Rcpp;
RcppExport SEXP gpuQR_magma(SEXP X_)
// Input
NumericMatrix X(X_);
// Initialize magma and cublas
magma_init();
cublasInit();
// Declare variables
int info, lwork, n_rows = X.nrow(), n_cols = X.ncol(), min_mn = min(n_rows, n_cols);
double tmp[1];
NumericVector scale(min_mn);
// Query workspace size
magma_dgeqrf(n_rows, n_cols, &(X[0]), n_rows, &(scale[0]), &(work[0]), -1, &info);
lwork = work[0];
NumericVector work(lwork);
// Run QR decomposition
magma_dgeqrf(n_rows, n_cols, &(X[0]), n_rows, &(scale[0]), &(work[0]), lwork, &info);
// Scale factor result
for(int ii = 1; ii < n_rows; ii++)
for(int jj = 0; jj < n_cols; jj++)
if(ii > jj) X[ii + jj * n_rows] *= scale[jj];
// Shutdown magma and cublas
magma_finalize();
cublasShutdown();
// Output
return wrap(X);
可以使用以下命令将文件从 R 编译到共享库中:
library(Rcpp)
PKG_LIBS <- sprintf('-Wl,-rpath,/usr/local/magma/lib -L/usr/local/magma/lib -lmagma /usr/local/magma/lib/libmagma.a -Wl,-rpath,/usr/local/cuda-5.5/lib64 %s', Rcpp:::RcppLdFlags())
PKG_CPPFLAGS <- sprintf('-DADD_ -DHAVE_CUBLAS -I/usr/local/magma/include -I/usr/local/cuda-5.5/include %s', Rcpp:::RcppCxxFlags())
Sys.setenv(PKG_LIBS = PKG_LIBS , PKG_CPPFLAGS = PKG_CPPFLAGS)
R <- file.path(R.home(component = 'bin'), 'R')
file <- '/path/gpuQR_magma.cpp'
cmd <- sprintf('%s CMD SHLIB %s', R, paste(file, collapse = ' '))
system(cmd)
现在可以在 R 中调用共享库。将结果与 R 的 qr() 进行比较得出:
dyn.load('/path/gpuQR_magma.so')
set.seed(100)
n_row <- 3; n_col <- 3
A <- matrix(rnorm(n_row * n_col), n_row, n_col)
qr(A)$qr
[,1] [,2] [,3]
[1,] 0.5250957 -0.8666925 0.8594266
[2,] -0.2504899 -0.3878643 -0.1277838
[3,] 0.1502909 0.4742033 -0.8804247
.Call('gpuQR_magma', A)
[,1] [,2] [,3]
[1,] 0.5250957 -0.8666925 0.8594266
[2,] -0.2504899 -0.3878643 -0.1277838
[3,] 0.1502909 0.4742033 -0.8804247
以下是使用具有 960 个 CUDA 内核和 OpenBLAS 的 NVIDIA GeForce GTX 675MX GPU 的基准测试结果:
n_row <- 3000; n_col <- 3000
A <- matrix(rnorm(n_row * n_col), n_row, n_col)
B <- A; dim(B) <- NULL
res <- benchmark(.Call('gpuQR_magma', A), .Call('gpuQR_cula', B, n_row, n_col), qr(A), columns = c('test', 'replications', 'elapsed', 'relative'), order = 'relative')
test replications elapsed relative
2 .Call("gpuQR_cula", B, n_row, n_col) 100 18.704 1.000
1 .Call("gpuQR_magma", A) 100 70.461 3.767
3 qr(A) 100 985.411 52.685
似乎 MAGMA 比 culatools 慢一点(在这个例子中)。但是,MAGMA 带有内存不足的实现,鉴于只有 1Gb 的 GPU 内存,我非常重视这一点。
【讨论】:
Rcpp 作者鼓励有像你这样的好例子的人在gallery.rcpp.org 上发帖,分享像你这样的好例子。从我的观点来看,你的东西应该在那里占有一席之地......让我们看看他们的想法 现有的 CRAN 包如 WideLm 使用 Rcpp 从 R 获取基于 C 的 CUDA API;可以轻松地为 Magma 做同样的事情。 实际上有magma R-package。但它不使用 Rcpp 。我发布了我的代码以供以后参考。 MAGMA 的文档在与 Rcpp 一起使用时非常薄,甚至更薄。 多哈。应该记得 CRAN 包岩浆。 是的,我把那个包作为我的代码的起点。向 Rcpp 库提交 Rcpp/MAGMA 示例是否有意义?以上是关于MAGMA 和 Rcpp 用于 R 中的线性代数的主要内容,如果未能解决你的问题,请参考以下文章