[Rcpp-devel] Dispatch on sparse and dense matrices...

Kevin Ushey kevinushey at gmail.com
Sun Dec 8 23:17:49 CET 2013


I think the safest bet would be to handle dispatch at the R level, which
would then call your internal C++ implementations. At least, it would be
easier / cleaner. Otherwise, you're right, you need to perform an 'ad-hoc'
dispatch based on the internal type / class and such.

-Kevin


On Sun, Dec 8, 2013 at 2:08 PM, Søren Højsgaard <sorenh at math.aau.dk> wrote:

> Thanks, Dirk!
>
> Follow up: That also means (I guess) that to be on the safe side I should
> write a check-function that will check that the INTSXP and REALSXP really
> is a matrix, i.e. that it has a dim attribute; right? Or are there
> facilities for that?
>
> Cheers
> Søren
>
>
> -----Original Message-----
> From: Dirk Eddelbuettel [mailto:edd at debian.org]
> Sent: 8. december 2013 22:59
> To: Søren Højsgaard
> Cc: rcpp-devel at lists.r-forge.r-project.org
> Subject: Re: [Rcpp-devel] Dispatch on sparse and dense matrices...
>
>
> Hi Søren,
>
> On 8 December 2013 at 21:47, Søren Højsgaard wrote:
> | Dear all,
> |
> | I want to do some computation on a matrix, irrespective of whether it is
> a numeric matrix, an integer matrix or a sparse matrix (a dgCMatrix which
> can be handled with RcppEigen). For simplicity, I want to compute the sum
> of the elements.
> |
> | To do so I use a template
> |
> | template <typename TT>
> | SEXP do_compute( SEXP XX_ ){
> |   const TT X(as<TT>(XX_));
> |   return wrap(X.sum()); // some computation on X...
> | };
> |
> | I then dispatch depending on input type with
> |
> | // [[Rcpp::export]]
> | SEXP compute ( SEXP XX_ ){
> |   int type = TYPEOF(XX_) ;
> |   Rf_PrintValue(wrap(type));
> |   switch( type ){
> |   case 13 : return do_compute<MapMati>(XX_); // matrix - integer
> |   case 14 : return do_compute<MapMatd>(XX_); // matrix - double
> |   case 25 : return do_compute<MSpMat>(XX_); // dgCMatrix
> |   }
> |   return R_NilValue ;
> | }
> |
> | and the numbers 13, 14, 25 are disclosed to by simply printing the type.
> I know that for integers and doubles I can do something "nicer":
> |
> | // [[Rcpp::export]]
> | SEXP compute2 ( SEXP XX_ ){
> |   int type = TYPEOF(XX_) ;
> |   Rf_PrintValue(wrap(type));
> |   switch( type ){
> |   case INTSXP : return do_compute<MapMati>(XX_); // matrix - integer
> |   case REALSXP : return do_compute<MapMatd>(XX_); // matrix - double
> |   case 25 : return do_compute<MSpMat>(XX_); // dgCMatrix
> |   }
> |   return R_NilValue ;
> | }
> |
> | Questions:
> |
> | Is there a similar keyword for sparse matrices??
>
> I don't think so as sparse matrices do not have a first-class
> representation in a SEXP. So no TYPEOF(...) that uniquely identifies them.
>
> | If not, is the "code" 25 safe to use in an R-package??
>
> No, it simply says 'S4'. You can get that for lots of things which are not
> Sparse matrices. From Rinternals:
>
>   #define S4SXP       25    /* S4, non-vector */
>
> So you could replace '25' with 'S4SXP' but that does not protect you from
> S4 types which are not sparse matrices.
>
> | Are there more elegant ways of handling the dispatch ??
>
> Not sure. The TYPEOF() idiom used to be the standard; I think a few nicer
> C++ wrappers for is<> were added more recently.  Those will still
> C++ dispatch on
> base types, and sparse matrices are not base types.  So you may have to go
> via S4SXP _and_ checking that the required slots exist.  Which you could
> then encode is a isSpMat() helper ...
>
> Dirk
>
> | The source file is attached. Thanks in advance.
> |
> | Søren
> |
> |
> |
> | ----------------------------------------------------------------------
> | #include <RcppEigen.h>
> | //[[Rcpp::depends(RcppEigen)]]
> |
> | using namespace Rcpp;
> |
> | typedef Eigen::Map<Eigen::MatrixXd> MapMatd; typedef
> | Eigen::Map<Eigen::MatrixXi> MapMati; typedef
> | Eigen::MappedSparseMatrix<double> MSpMat;
> |
> | // [[Rcpp::export]]
> | SEXP getType ( SEXP XX_ ){
> |   int type = TYPEOF(XX_) ;
> |   return wrap(type);
> | }
> |
> | template <typename TT>
> | SEXP do_compute( SEXP XX_ ){
> |   const TT X(as<TT>(XX_));
> |   return wrap(X.sum()); // some computation on X...
> | };
> |
> | // [[Rcpp::export]]
> | SEXP compute ( SEXP XX_ ){
> |   int type = TYPEOF(XX_) ;
> |   Rf_PrintValue(wrap(type));
> |   switch( type ){
> |   case 13 : return do_compute<MapMati>(XX_); // matrix - integer
> |   case 14 : return do_compute<MapMatd>(XX_); // matrix - double
> |   case 25 : return do_compute<MSpMat>(XX_); // dgCMatrix
> |   }
> |   return R_NilValue ;
> | }
> |
> | // [[Rcpp::export]]
> | SEXP compute2 ( SEXP XX_ ){
> |   int type = TYPEOF(XX_) ;
> |   Rf_PrintValue(wrap(type));
> |   switch( type ){
> |   case INTSXP : return do_compute<MapMati>(XX_); // matrix - integer
> |   case REALSXP : return do_compute<MapMatd>(XX_); // matrix - double
> |   case 25 : return do_compute<MSpMat>(XX_); // dgCMatrix
> |   }
> |   return R_NilValue ;
> | }
> |
> | /*** R
> |
> | library(Matrix)
> | m <- matrix(1:9, nrow=3)
> | M <- as(m, "dgCMatrix")
> |
> | getType( m ); getType( M ); getType( 1 ); getType( 1L )
> |
> | compute( m )   # integer
> | compute( 1*m ) # double
> | compute( M )   # sparse dgCMatrix
> | */
> |
> |
> |
> | ----------------------------------------------------------------------
> | _______________________________________________
> | Rcpp-devel mailing list
> | Rcpp-devel at lists.r-forge.r-project.org
> | https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-deve
> | l
> --
> Dirk Eddelbuettel | edd at debian.org | http://dirk.eddelbuettel.com
> _______________________________________________
> Rcpp-devel mailing list
> Rcpp-devel at lists.r-forge.r-project.org
> https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.r-forge.r-project.org/pipermail/rcpp-devel/attachments/20131208/594996b1/attachment.html>


More information about the Rcpp-devel mailing list