MAGMA and Rcpp for linear algebra in R

I was wondering if anybody has tried to use Rcpp and MAGMA to accelerate linear algebra operations in R by using the CPU and GPU? I tried culatools last month and it worked with Rcpp (link), but culatools is a commercial product that costs money to get access to all functions.

Answers


It was pretty straightforward to use Rcpp and MAGMA after having tinkered around with culatools. Here is the .cpp file:

#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);
}

The file can be compiled from R into a shared library using:

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)

The shared library can now be called in R. Comparing the results with with R's qr() gives:

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

Below are the results from a benchmark using a NVIDIA GeForce GTX 675MX GPU with 960 CUDA cores and 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

Seems like as if MAGMA is a bit slower compared to culatools (in this example). However, MAGMA comes with out-of-memory implementations and that is something I value a lot given only 1Gb of GPU memory.


Need Your Help

Is there any indication when Android kills a service?

android service kill-process

I have a service running in its own process. It appears fine in the evening, but after I've gone to sleep I think Android takes its axe to it.

About UNIX Resources Network

Original, collect and organize Developers related documents, information and materials, contains jQuery, Html, CSS, MySQL, .NET, ASP.NET, SQL, objective-c, iPhone, Ruby on Rails, C, SQL Server, Ruby, Arrays, Regex, ASP.NET MVC, WPF, XML, Ajax, DataBase, and so on.