Estou desenvolvendo algumas funções emRcpp que operam embig.matrix objetos dobigmemory pacote. Esses objetos são passados paraRcpp ComoSEXP objetos que eu tenho que converter para umXPtr<BigMatrix>e depois para umMatrixAccessor objeto para acessar elementos da matriz.

Por exemplo, se eu quiser implementar uma função que é da 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;
}

Esta função funciona perfeitamente, desde que obig.matrix objeto em R é preenchido com duplas.

No entanto, se você chamar essa função em uma matriz inteira (por exemplo,diag(as.big.matrix(matrix(1:9, 3))@address)), Você recebe lixo como resultado, porque oMatrixAccessor foi inicializado como<double>.

Internamente,big.matrix objetos podem vir em quatro tipos:

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
}

Como tudo o que estamos fazendo é acessar os elementos da matriz, odiag A função deve poder lidar com cada um desses tipos. Mas, por enquanto, já que nossa assinatura de função éNumericVector, Ignorarei matrizes de caracteres.

Para lidar com isso, achei que poderia simplesmente lançar uma instrução switch, inicializando o correspondentemat com o tipo apropriado em tempo de execução:

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

No entanto, isso resulta em erros do compilador, porque estou redefinindomat depois de ter sido declarado já no primeirocase.

A única maneira de ver isso é escrever três diferentesdiag funções, uma para cada tipo, cujo código é o mesmo, com exceção da inicialização demat. Existe uma maneira melhor?