Using RViennaCL for matrix multiplication

Question

I'm a newbie using the RViennaCL package, trying to do GPU matrix multiplication.

Trying a simple case I created a cpp file inside the src folder within the package I created (gpuUtils), with this code:

#include "viennacl/ocl/backend.hpp"
#include "viennacl/linalg/prod.hpp"
#include "viennacl/tools/random.hpp"
#include "viennacl/matrix.hpp"


//' PU matrix multiplication
//'
//' @export
// [[Rcpp::export]]
SEXP matMult(){

  viennacl::tools::uniform_random_numbers<float> randomNumber;

  viennacl::matrix<float> vcl_A(400, 400);
  viennacl::matrix<float> vcl_B(400, 400);

  for (unsigned int i = 0; i < vcl_A.size1(); ++i)
    for (unsigned int j = 0; j < vcl_A.size2(); ++j)
      vcl_A(i,j) = randomNumber();
  for (unsigned int i = 0; i < vcl_B.size1(); ++i)
    for (unsigned int j = 0; j < vcl_B.size2(); ++j)
      vcl_B(i,j) = randomNumber();
  viennacl::matrix<float> vcl_C = viennacl::linalg::prod(vcl_A, vcl_B);
  return Rcpp::wrap(vcl_C);
}

Then I build the package:

$ R CMD build gpuUtils
* checking for file ‘gpuUtils/DESCRIPTION’ ... OK
* preparing ‘gpuUtils’:
* checking DESCRIPTION meta-information ... OK
* cleaning src
* checking for LF line-endings in source and make files and shell scripts
* checking for empty or unneeded directories
* building ‘gpuUtils_0.0.0.9000.tar.gz’

But when trying to install it I get:

$ R CMD INSTALL gpuUtils_0.0.0.9000.tar.gz
* installing to library ‘/home/usr/R’
* installing *source* package ‘gpuUtils’ ...
** libs
g++  -I/opt/ohpc/pub/libs/gnu7/R/3.4.2/lib64/R/include -DNDEBUG -I/usr/local/cuda-9.0/include/ -I"/home/usr/R/Rcpp/include" -I"/home/usr/R/RViennaCL/include" -I/usr/local/include   -fpic  -g -O2  -c gpuUtils.cpp -o gpuUtils.o
In file included from /home/usr/R/Rcpp/include/RcppCommon.h:195:0,
                 from /home/usr/R/Rcpp/include/Rcpp.h:27,
                 from /home/usr/R/RViennaCL/include/viennacl/ocl/forwards.h:29,
                 from /home/usr/R/RViennaCL/include/viennacl/ocl/context.hpp:36,
                 from /home/usr/R/RViennaCL/include/viennacl/ocl/backend.hpp:26,
                 from gpuUtils.cpp:1:
/home/usr/R/Rcpp/include/Rcpp/internal/wrap.h: In instantiation of ‘SEXPREC* Rcpp::internal::wrap_dispatch_unknown_iterable(const T&, Rcpp::traits::false_type) [with T = viennacl::matrix<float, viennacl::row_major>; SEXP = SEXPREC*; Rcpp::traits::false_type = Rcpp::traits::integral_constant<bool, false>]’:
/home/usr/R/Rcpp/include/Rcpp/internal/wrap.h:732:43:   required from ‘SEXPREC* Rcpp::internal::wrap_dispatch_unknown(const T&, Rcpp::traits::false_type) [with T = viennacl::matrix<float, viennacl::row_major>; SEXP = SEXPREC*; Rcpp::traits::false_type = Rcpp::traits::integral_constant<bool, false>]’
/home/usr/R/Rcpp/include/Rcpp/internal/wrap.h:770:41:   required from ‘SEXPREC* Rcpp::internal::wrap_dispatch_eigen(const T&, Rcpp::traits::false_type) [with T = viennacl::matrix<float, viennacl::row_major>; SEXP = SEXPREC*; Rcpp::traits::false_type = Rcpp::traits::integral_constant<bool, false>]’
/home/usr/R/Rcpp/include/Rcpp/internal/wrap.h:787:39:   required from ‘SEXPREC* Rcpp::internal::wrap_dispatch_unknown_importable(const T&, Rcpp::traits::false_type) [with T = viennacl::matrix<float, viennacl::row_major>; SEXP = SEXPREC*; Rcpp::traits::false_type = Rcpp::traits::integral_constant<bool, false>]’
/home/usr/R/Rcpp/include/Rcpp/internal/wrap.h:807:52:   required from ‘SEXPREC* Rcpp::internal::wrap_dispatch(const T&, Rcpp::traits::wrap_type_unknown_tag) [with T = viennacl::matrix<float, viennacl::row_major>; SEXP = SEXPREC*]’
/home/usr/R/Rcpp/include/Rcpp/internal/wrap_end.h:30:38:   required from ‘SEXPREC* Rcpp::wrap(const T&) [with T = viennacl::matrix<float, viennacl::row_major>; SEXP = SEXPREC*]’
gpuUtils.cpp:25:26:   required from here
/home/usr/R/Rcpp/include/Rcpp/internal/wrap.h:523:17: error: static assertion failed: cannot convert type to SEXP
                 static_assert(!sizeof(T), "cannot convert type to SEXP");
                 ^~~~~~~~~~~~~
make: *** [gpuUtils.o] Error 1
ERROR: compilation failed for package ‘gpuUtils’
* removing ‘/home/usr/R/gpuUtils’
* restoring previous ‘/home/usr/R/gpuUtils’

Which indicates that viennacl::matrix need its own wrap function to be written.

So I tried to follow this very nice vignette for custom templating as and wrap functions within Rcpp, adapting it to the viennacl::matrix case.

Here's my cpp code:

#include <RcppCommon.h>
#include "viennacl/ocl/backend.hpp"
#include "viennacl/matrix.hpp"


// [[Rcpp::plugins("cpp11")]]
// Provide Forward Declarations
namespace Rcpp {
  namespace traits{
    // Setup non-intrusive extension via template specialization for
    // 'viennacl' class

    // Support for wrap
    template <typename T> SEXP wrap(const viennacl::matrix<T> & obj);

    // Support for as<T>
    template <typename T> class Exporter< viennacl::matrix<T> >;
  }
}

#include <Rcpp.h>

// Define template specializations for as<> and wrap
namespace Rcpp {

  namespace traits{

  // Defined wrap case
  template <typename T> SEXP wrap(const viennacl::matrix<T> & obj){
    const int RTYPE = Rcpp::traits::r_sexptype_traits<T>::rtype ;

    return Rcpp::Matrix< RTYPE >(obj.begin(), obj.end());
  };


  // Defined as< > case
  template<typename T> class Exporter< viennacl::matrix<T> > {
    typedef typename viennacl::matrix<T> OUT ;

    // Convert the type to a valid rtype.
    const static int RTYPE = Rcpp::traits::r_sexptype_traits< T >::rtype ;
    Rcpp::Matrix<RTYPE> mat;

    public:
      Exporter(SEXP x) : mat(x) {
        if (TYPEOF(x) != RTYPE)
          throw std::invalid_argument("Wrong R type for mapped 2D array");
      }
      OUT get() {

       // Need to figure out a way to perhaps do a pointer pass?
        OUT x(mat.size());

        std::copy(mat.begin(), mat.end(), mat.begin());

        return x;
      }
    };
  }
}

//' viennacl matrix templating
//'
//' @export
// [[Rcpp::export]]
void matTemplate(Rcpp::NumericMatrix in_mat){
  viennacl::matrix<double> viennacl_mat = Rcpp::as< viennacl::matrix<double> >(in_mat);
  Rcpp::NumericMatrix out_mat = Rcpp::wrap(viennacl_mat);
}

Builds fine:

$ R CMD build gpuUtils
* checking for file ‘gpuUtils/DESCRIPTION’ ... OK
* preparing ‘gpuUtils’:
* checking DESCRIPTION meta-information ... OK
* cleaning src
* checking for LF line-endings in source and make files and shell scripts
* checking for empty or unneeded directories
* building ‘gpuUtils_0.0.0.9000.tar.gz’

But still getting the same error:

$ R CMD INSTALL gpuUtils_0.0.0.9000.tar.gz
* installing to library ‘/home/usr/R’
* installing *source* package ‘gpuUtils’ ...
** libs
g++  -I/opt/ohpc/pub/libs/gnu7/R/3.4.2/lib64/R/include -DNDEBUG -I/usr/local/cuda-9.0/include/ -I"/home/usr/R/Rcpp/include" -I"/home/usr/R/RViennaCL/include" -I/usr/local/include   -fpic  -g -O2  -c gpuUtils.cpp -o gpuUtils.o
In file included from /home/usr/R/Rcpp/include/RcppCommon.h:195:0,
                 from gpuUtils.cpp:1:
/home/usr/R/Rcpp/include/Rcpp/internal/wrap.h: In instantiation of ‘SEXPREC* Rcpp::internal::wrap_dispatch_unknown_iterable(const T&, Rcpp::traits::false_type) [with T = viennacl::matrix<double, viennacl::row_major>; SEXP = SEXPREC*; Rcpp::traits::false_type = Rcpp::traits::integral_constant<bool, false>]’:
/home/usr/R/Rcpp/include/Rcpp/internal/wrap.h:732:43:   required from ‘SEXPREC* Rcpp::internal::wrap_dispatch_unknown(const T&, Rcpp::traits::false_type) [with T = viennacl::matrix<double, viennacl::row_major>; SEXP = SEXPREC*; Rcpp::traits::false_type = Rcpp::traits::integral_constant<bool, false>]’
/home/usr/R/Rcpp/include/Rcpp/internal/wrap.h:770:41:   required from ‘SEXPREC* Rcpp::internal::wrap_dispatch_eigen(const T&, Rcpp::traits::false_type) [with T = viennacl::matrix<double, viennacl::row_major>; SEXP = SEXPREC*; Rcpp::traits::false_type = Rcpp::traits::integral_constant<bool, false>]’
/home/usr/R/Rcpp/include/Rcpp/internal/wrap.h:787:39:   required from ‘SEXPREC* Rcpp::internal::wrap_dispatch_unknown_importable(const T&, Rcpp::traits::false_type) [with T = viennacl::matrix<double, viennacl::row_major>; SEXP = SEXPREC*; Rcpp::traits::false_type = Rcpp::traits::integral_constant<bool, false>]’
/home/usr/R/Rcpp/include/Rcpp/internal/wrap.h:807:52:   required from ‘SEXPREC* Rcpp::internal::wrap_dispatch(const T&, Rcpp::traits::wrap_type_unknown_tag) [with T = viennacl::matrix<double, viennacl::row_major>; SEXP = SEXPREC*]’
/home/usr/R/Rcpp/include/Rcpp/internal/wrap_end.h:30:38:   required from ‘SEXPREC* Rcpp::wrap(const T&) [with T = viennacl::matrix<double, viennacl::row_major>; SEXP = SEXPREC*]’
gpuUtils.cpp:68:56:   required from here
/home/usr/R/Rcpp/include/Rcpp/internal/wrap.h:523:17: error: static assertion failed: cannot convert type to SEXP
                 static_assert(!sizeof(T), "cannot convert type to SEXP");
                 ^~~~~~~~~~~~~
make: *** [gpuUtils.o] Error 1
ERROR: compilation failed for package ‘gpuUtils’
* removing ‘/home/usr/R/gpuUtils’
* restoring previous ‘/home/usr/R/gpuUtils’

Answer 1

Short answer: yes. RViennaCl only provides the ViennaCl headers but no as() and wrap() functions for communication between R and C++ data structures. And the gpuR package, which uses RViennaCl, uses a different approach by exposing suitable objects to R with appropriate functions acting on these objects. Of course, you can use these objects to do matrix multiplication on the GPU.

So you would have to write these functions yourself if you are interested in communication between R and ViennaCl data structures. This is not a simple task, though, since RViennaCL patches the ViennaCL headers to include Rcpp.h (e.g. https://github.com/cdeterman/RViennaCL/blob/master/inst/include/viennacl/ocl/backend.hpp#L29). This makes it impossible to use the usual extension mechanisms as you are trying above. I think it would be sufficient for RViennaCL to include RcppCommon.h, this way making extensions possible. But I have not tried that. I have filed https://github.com/cdeterman/RViennaCL/issues/4 for this.

Instead of writing as() and wrap() for ViennaCL you could also combine available converters with a package that converts to and from R data structures. So in a package that uses LinkingTo: Rcpp, RViennaCL, RcppEigen this function works for me:

#define VIENNACL_HAVE_EIGEN
#include <RcppEigen.h>
#include <viennacl/ocl/backend.hpp>
#include <viennacl/matrix.hpp>


//' viennacl matrix crossprod
//'
//' @export
// [[Rcpp::export]]
Eigen::MatrixXd vclCrossprod(const Eigen::MatrixXd &in_densematrix){
  int rows = in_densematrix.rows();
  int cols = in_densematrix.cols();

  viennacl::matrix<double>  viennacl_densematrix(rows, cols);
  viennacl::copy(in_densematrix,  viennacl_densematrix);
  viennacl::matrix<double> viennacl_result = viennacl::linalg::prod(
    viennacl::trans(viennacl_densematrix),
    viennacl_densematrix);

  Eigen::MatrixXd eigen_result(cols, cols);
  viennacl::copy(viennacl_result,  eigen_result);
  return eigen_result;
}

However, for the examples I tried, this approach is not really efficient. This might change if you are doing something more complicated than matrix multiplication.

Alternatively you could try RcppArrayFire, which I am working on, if you are mainly interested in GPGPU and not ViennaCl in particular.