An example using RViennaCL for matrix multiplication

[ndrubins]

Hi @cdeterman ,

Would it be possible to get an example code that does matrix multiplication using RViennaCL?

Thanks a lot,

[cdeterman]

Hi @ndrubins, currently the best I have available would be in my gpuR package. The R matrix needs to first be copied to a viennacl::matrix object. It is my intent to make some nice copy function for this but it does not exist yet. I was initially leveraging the NumericMatrix class from Rcpp but that would only work for the double type so I will need to write with the base R logic. Regardless, if you want to do it, you will likely need two steps.

1. Copy the matrix to a Eigen::Matrix.
2. Copy the Eigen::Matrix to a viennacl::matrix

These two steps are demonstrated in this section of gpuR here

3. Perform you matrix multiplication The template I use in gpuR can be found here.

4. Make sure to copy the data back to the host (see last line of step 3) whenever you want it returned back to R. Unless you want it to remain on the device.

Naturally, all of this assumes you are trying to leverage a GPU. If you are using the OpenMP backend, you shouldn't need to worry about copying back and forth with host->device->host. But I am not as familiar with that and you would need to consult the ViennaCL documentation.

I hope this helps. I want to work on RViennaCL and make it more usable beyond the wonderful headers but I only have so much time and other things take priority currently.

[cdeterman]

@ndrubins I would like to add though if you make an attempt at this and run in to problems. Feel free to upload your code to git and I will help when I can.

[ndrubins]

Thanks a lot.

I started by trying to write a simple matrix multiplication function in a cpp file (which is in a package I called gpuUtils):

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


//' GPU 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);
}

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 seems like a viennacl::matrix cannot be simply wrapped by Rcpp::wrap:

$ R CMD INSTALL gpuUtils_0.0.0.9000.tar.gz
* installing to library ‘/home/usr/R’
* installing *source* package ‘gpuUtils’ ...
** libs
g++  -I/sw/R/R-3.4.3-install/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’

So I tried to define as and wrap templates for viennacl::matrix (following this vignette), and this time simply trying a function that converts an Rcpp::NumericMatrix to a viennacl::matrix<double>:

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

//' convert and Rcpp::NumericMatrix to a viennacl::matrix
//'
//' @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);
}

But still failing:

$ R CMD INSTALL gpuUtils_0.0.0.9000.tar.gz
* installing to library ‘/home/usr/R’
* installing *source* package ‘gpuUtils’ ...
** libs
g++  -I/sw/R/R-3.4.3-install/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’