MAGMA и Rcpp для линейной алгебры в R

Question

Мне было интересно, пытался ли кто-нибудь использовать Rcpp и MAGMA для ускорения операций линейной алгебры в R с использованием процессора и графического процессора? Я пробовал culatools в прошлом месяце, и он работал с Rcpp (link), но culatools - это коммерческий продукт, который стоит денег, чтобы получить доступ ко всем функциям.

Answer 1

Это было довольно проста в использовании Rcpp и МАГМА после того, как возился вокруг с culatools. Вот файл .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 «s 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 

Ниже приведены результаты из теста с помощью NVIDIA GeForce GTX 675MX GPU с 960 CUDA ядрами и OpenBLAS:

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 поставляется с реализациями вне памяти, и это то, что я очень ценю, учитывая только 1 Гб памяти GPU.