Ich entwickle einige Funktionen inRcpp die weiterarbeitenbig.matrix Objekte aus dembigmemory Paket. Diese Objekte werden an übergebenRcpp wieSEXP Objekte, die ich dann zu einemXPtr<BigMatrix>und dann zu aMatrixAccessor Objekt, um auf Elemente der Matrix zuzugreifen.

Wenn ich zum Beispiel eine Funktion implementieren möchte, die die Diagonale erhält:

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

Diese Funktion funktioniert einwandfrei, sofern diebig.matrix Objekt in R ist mit Doppeln gefüllt.

Wenn Sie diese Funktion jedoch für eine Ganzzahlmatrix aufrufen (z.diag(as.big.matrix(matrix(1:9, 3))@address)), Sie bekommen als Ergebnis Müll, weil dieMatrixAccessor wurde initialisiert als<double>.

Im Inneren,big.matrix Es gibt vier Arten von Objekten:

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
}

Da wir nur auf die Elemente der Matrix zugreifen, wird diediag Funktion sollte in der Lage sein, mit jedem dieser Typen umzugehen. Aber vorerst ist da unsere FunktionsunterschriftNumericVectorIch werde Zeichenmatrizen ignorieren.

Um das zu handhaben, dachte ich, ich könnte einfach eine switch-Anweisung einbauen und die entsprechende initialisierenmat mit dem entsprechenden Typ zur Laufzeit:

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

Dies führt jedoch zu Compilerfehlern, da ich sie neu definieremat nachdem es schon im ersten deklariert wurdecase.

Die einzige Möglichkeit, dies zu tun, besteht darin, drei verschiedene zu schreibendiag Funktionen, eine für jeden Typ, deren Code mit Ausnahme der Initialisierung von gleich istmat. Gibt es einen besseren Weg?