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

Question

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

Answer 1

Использовать 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

Ниже приведены результаты теста производительности графического процессора NVIDIA GeForce GTX 675MX с ядрами 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 ГБ памяти графического процессора.