RcppArmadillo中稀疏和密集矩阵的模板函数

Question

我试图定义一个模板函数，它可以使用RcppArmadillo处理稀疏和密集矩阵输入。我发送一个密集或稀疏矩阵C++和回R的非常简单的情况下，像这样的工作：

library(inline); library(Rcpp); library(RcppArmadillo) 

sourceCpp(code = " 
// [[Rcpp::depends(RcppArmadillo)]] 
#include <RcppArmadillo.h> 
using namespace Rcpp ; 
using namespace arma ; 

template <typename T> T importexport_template(const T X) { 
    T ret = X ; 
    return ret ; 
}; 

//[[Rcpp::export]] 
SEXP importexport(SEXP X) { 
    return wrap(importexport_template(X)) ; 
}") 

library(Matrix) 
X <- diag(3) 
X_sp <- as(X, "dgCMatrix") 

importexport(X) 
##  [,1] [,2] [,3] 
##[1,] 1 0 0 
##[2,] 0 1 0 
##[3,] 0 0 1 
importexport(X_sp) 
##3 x 3 sparse Matrix of class "dgCMatrix" 
##   
##[1,] 1 . . 
##[2,] . 1 . 
##[3,] . . 1 

，我理解这意味着模板基本工作原理（即密集R-矩阵通过对Rcpp::as的隐式调用，以及相应的暗指Rcpp:wrap然后做正确的事情并返回稠密以用于稀疏稀疏，稀疏R矩阵变为arma::sp_mat-对象。

的实际功能我尝试写需求，当然多个参数，而这也正是我失败 - 做这样的事情：

sourceCpp(code =" 
// [[Rcpp::depends(RcppArmadillo)]] 
#include <RcppArmadillo.h> 

using namespace Rcpp ; 
using namespace arma ; 

template <typename T> T scalarmult_template(const T X, double scale) { 
    T ret = X * scale; 
    return ret; 
}; 

//[[Rcpp::export]] 
SEXP scalarmult(SEXP X, double scale) { 
    return wrap(scalarmult_template(X, scale)) ; 
}") 

失败，因为编译器不知道如何在编译时解决*为SEXPREC* const。 所以我想我需要的东西，如开关语句in this Rcpp Gallery snippet正确分派到特定的模板功能，但我不知道该怎么写了，似乎比INTSXP等

比较复杂，我想我知道如何种访问类型我需要这样一个switch语句，例如：

sourceCpp(code =" 
// [[Rcpp::depends(RcppArmadillo)]] 
#include <RcppArmadillo.h> 

using namespace Rcpp ; 
using namespace arma ; 

//[[Rcpp::export]] 
SEXP printtype(SEXP Xr) { 
    Rcpp::Rcout << TYPEOF(Xr) << std::endl ; 
    return R_NilValue; 
}") 
printtype(X) 
##14 
##NULL 
printtype(X_sp) 
##25 
##NULL 

但我不明白如何从那里继续。 scalarmult_template适用于稀疏矩阵和密集矩阵的版本是什么样的？

Answer 1

基于@ KevinUshey的评论回答我自己的问题。我做矩阵乘法为三种情况：密密，疏密，以及“indMatrix” -dense：

library(inline) 
library(Rcpp) 
library(RcppArmadillo) 
library(Matrix) 
library(rbenchmark) 

sourceCpp(code=" 
// [[Rcpp::depends(RcppArmadillo)]] 
#include <RcppArmadillo.h> 

using namespace Rcpp ; 
using namespace arma ; 

arma::mat matmult_sp(const arma::sp_mat X, const arma::mat Y){ 
    arma::mat ret = X * Y; 
    return ret; 
}; 
arma::mat matmult_dense(const arma::mat X, const arma::mat Y){ 
    arma::mat ret = X * Y; 
    return ret; 
}; 
arma::mat matmult_ind(const SEXP Xr, const arma::mat Y){ 
    // pre-multplication with index matrix is a permutation of Y's rows: 
    S4 X(Xr); 
    arma::uvec perm = X.slot("perm"); 
    arma::mat ret = Y.rows(perm - 1); 
    return ret; 
}; 

//[[Rcpp::export]] 
arma::mat matmult_cpp(SEXP Xr, const arma::mat Y) { 
    if (Rf_isS4(Xr)) { 
     if(Rf_inherits(Xr, "dgCMatrix")) { 
      return matmult_sp(as<arma::sp_mat>(Xr), Y) ; 
     } ; 
     if(Rf_inherits(Xr, "indMatrix")) { 
      return matmult_ind(Xr, Y) ; 
     } ; 
     stop("unknown class of Xr") ; 
    } else { 
     return matmult_dense(as<arma::mat>(Xr), Y) ; 
    } 
}") 

n <- 10000 
d <- 20 
p <- 30 

X <- matrix(rnorm(n*d), n, d) 
X_sp <- as(diag(n)[,1:d], "dgCMatrix") 
X_ind <- as(sample(1:d, n, rep=TRUE), "indMatrix") 
Y <- matrix(1:(d*p), d, p) 

matmult_cpp(as(X_ind, "ngTMatrix"), Y) 
## Error: unknown class of Xr 

all.equal(X%*%Y, matmult_cpp(X, Y)) 
## [1] TRUE 

all.equal(as.vector(X_sp%*%Y), 
      as.vector(matmult_cpp(X_sp, Y))) 
## [1] TRUE 

all.equal(X_ind%*%Y, matmult_cpp(X_ind, Y)) 
## [1] TRUE 

---

编辑：这已经变成一个Rcpp Gallery post。