Initialize a variable with a different type based on the switch statement

Question

I am developing some functions in Rcpp that work with big.matrix objects from the bigmemory package. These objects are passed to Rcpp as SEXP objects, which then need to be passed to XPtr<BigMatrix> , and then to the MatrixAccessor object to access the matrix elements.

For example, if I want to implement a function that receives a diagonal:

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

This function works beautifully if the big.matrix object in R is filled with doubles.

However, if you call this function on an integer matrix (for example, diag(as.big.matrix(matrix(1:9, 3))@address) ), you get garbage as a result because MatrixAccessor was initialized as <double> .

Internally, big.matrix objects can have four types:

 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 } 

Since all we do is access the matrix elements, the diag function should be able to handle each of these types. But for now, since our signature NumericVector , I will ignore the character matrices.

To handle this, I decided that I could just specify the switch statement by initializing the corresponding mat with the appropriate type at runtime:

 // [[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; } 

However, this leads to compiler errors, because I override mat after it was declared already in the first case .

The only way I can do this is to write three different diag functions, one for each type, whose code is the same, except for the initialization mat . Is there a better way?

Answer 1

Thanks to Kevin Ears' pointers to the separation of control room and function logic. The required template code is sufficiently different for Kevin's suggestions to guarantee his own answer.

To write a function that usually works for all big.matrix types, the template looks like this:

 #include <Rcpp.h> using namespace Rcpp; // [[Rcpp::depends(BH, bigmemory)]] #include <bigmemory/MatrixAccessor.hpp> // Generic Implementation of GetDiag: template <typename T> NumericVector GetDiag_impl(XPtr<BigMatrix> xpMat, MatrixAccessor<T> mat) { NumericVector diag(xpMat->nrow()); // Assume matrix is square for (unsigned int i = 0; i < xpMat->nrow(); i++) { diag[i] = mat[i][i]; } return diag; } // Dispatch code // [[Rcpp::export]] NumericVector GetDiag(SEXP pMat) { XPtr<BigMatrix> xpMat(pMat); switch(xpMat->matrix_type()) { case 1: return GetDiag_impl(xpMat, MatrixAccessor<char>(*xpMat)); case 2: return GetDiag_impl(xpMat, MatrixAccessor<short>(*xpMat)); case 4: return GetDiag_impl(xpMat, MatrixAccessor<int>(*xpMat)); case 8: return GetDiag_impl(xpMat, MatrixAccessor<double>(*xpMat)); default: // This should be impossible to reach, but shuts up the compiler throw Rcpp::exception("Unknown type of big.matrix detected! Aborting."); } } 

You do not need to create a template of the GetDiag_impl return type: all four big.matrix types big.matrix stored as numeric in R (see this answer for a discussion on the 'char' big.matrix objects).

Answer 2

In these cases, you usually want to separate the logic: first, you have a submit function that sorts the SEXP type into some type of compilation time and a separate (template) function that handles the actual work. Some (incomplete) code example:

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