rcpp - 根据switch语句初始化具有不同类型的变量

Question

我正在Rcpp中开发一些对big.matrix包中的bigmemory对象进行操作的函数。这些对象作为Rcpp对象传递给SEXP，然后我必须将它们转换为XPtr<BigMatrix>，然后转换为MatrixAccessor对象以访问矩阵的元素。

例如，如果我想实现一个以对角线表示的函数:

#include <Rcpp.h>
using namespace Rcpp;

// [[Rcpp::depends(BH, bigmemory)
#include <bigmemory/MatrixAccessor.hpp>

#include <numeric>

// [[Rcpp::export]]
NumericVector GetDiag(SEXP pmat) {
  XPtr<BigMatrix> xpMat(pMat); // allows you to access attributes
  MatrixAccessor<double> mat(*xpMat); // allows you to access matrix elements

  NumericVector diag(xpMat->nrow()); // Assume matrix is square
  for (int i = 0; i < xpMat->nrow(); i++) { 
    diag[i] = mat[i][i];
  }
  return diag;
}

只要R中的 big.matrix对象用 double 填充，此函数就可以很好地工作。

但是，如果您在整数矩阵(例如 diag(as.big.matrix(matrix(1:9, 3))@address))上调用此函数，则结果会产生垃圾，因为 MatrixAccessor已初始化为 <double>。

在内部， big.matrix对象可以分为四种类型:

void typeof(SEXP pMat) {
  XPtr<BigMatrix> xpMat(pMat);
  int type = xpMat->matrix_type();
  type == 1 // char
  type == 2 // short
  type == 4 // int
  type == 8 // double
}

由于我们要做的只是访问矩阵的元素，因此 diag函数应该能够处理所有这些类型。但是现在，由于我们的函数签名是 NumericVector，所以我将忽略字符矩阵。

为了解决这个问题，我认为我可以抛出一个switch语句，在运行时用适当的类型初始化对应的 mat:

// [[Rcpp::export]]
NumericVector GetDiag(SEXP pmat) {
  XPtr<BigMatrix> xpMat(pMat);
  // Determine the typeof(pmat), and initialize accordingly:
  switch(xpMat->matrix_type()) {
    case == 1:
    {
       // Function expects to return a NumericVector.
       throw; 
    }
    case == 2:
    {
      MatrixAccessor<short> mat(*xpMat);
      break;
    }
    case == 4:
    {
      MatrixAccessor<int> mat(*xpMat);
      break;
    }
    case == 8:
    {
      MatrixAccessor<double> mat(*xpMat);
    }
  }
  MatrixAccessor<double> mat(*xpMat); // allows you to access matrix elements

  NumericVector diag(xpMat->nrow()); // Assume matrix is square
  for (int i = 0; i < xpMat->nrow(); i++) { 
    diag[i] = mat[i][i];
  }
  return diag;
}

但是，这会导致编译器错误，因为在第一个 mat中已经声明了 case之后，我要重新定义它。

我看到的唯一方法是编写三种不同的 diag函数，每种函数一个，除 mat的初始化外，其代码均相同。有没有更好的办法？

Answer 1

在这些情况下，通常需要分离逻辑:首先，您具有一个将SEXP类型编码为某种编译时类型的调度函数，以及一个用于处理实际工作的单独(模板化)函数。一些(不完整的)示例代码:

#include <Rcpp.h>
using namespace Rcpp;

// the actual generic implementation
template <typename T>
T GetDiag_impl(T const& pMat) {
  // logic goes here
}

// the dispatch code

// [[Rcpp::export]]
SEXP GetDiag(SEXP pMat) {
  switch (TYPEOF(pMat)) {
  case INTSXP: return GetDiag_impl<IntegerMatrix>(pMat);
  case REALSXP: return GetDiag_impl<NumericMatrix>(pMat);
  <...>
  }
}