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 矩阵变成了 arma::mat，而稀疏的 R 矩阵通过隐式调用变成了 arma::sp_mat-object Rcpp::as 和相应的隐含 Rcpp:wraps 然后也做正确的事情，并为密集返回密集，为稀疏返回稀疏）。

我尝试编写的实际函数当然需要多个参数，这就是我失败的地方 - 执行以下操作：

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 解析 *。 所以我想我需要类似 switch 语句 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 的评论回答我自己的问题。我对 3 种情况进行矩阵乘法：dense-dense、sparse-dense 和“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。