MAGMA et CRPP pour l'algèbre linéaire dans R

Question

Je me demandais si quelqu'un a essayé d'utiliser Rcpp et MAGMA pour accélérer les opérations d'algèbre linéaire dans R en utilisant le CPU et le GPU? J'ai essayé culatools le mois dernier et cela a fonctionné avec Rcpp (link), mais culatools est un produit commercial qui coûte de l'argent pour avoir accès à toutes les fonctions.

Answer 1

Il est assez simple à utiliser CRPP et MAGMA après avoir bricolé autour avec culatools. Voici le fichier .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); 
} 

Le fichier peut être compilé à partir de R dans une bibliothèque partagée en utilisant:

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) 

La bibliothèque partagée peut être appelée de R. En comparant les résultats obtenus avec des R de qr() donne:

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 

Voici les résultats d'une référence en utilisant un GPU NVIDIA GeForce GTX 675MX avec 960 cœurs CUDA et 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 

On dirait que si MAGMA est un peu plus lent par rapport à culatools (dans cet exemple). Cependant, MAGMA est livré avec des implémentations de mémoire insuffisante et c'est quelque chose que j'estime beaucoup avec seulement 1Go de mémoire GPU.