Casting a BigMatrix/array to a Armadillo matrix

Question

I have a big.matrix that I want to cast to an arma::Mat so that I can use the linear algebra functionality of Armadillo.

However, I can't seem to get the cast to work.

As far as I can gather from reading, both are internally stored in column major format, and the actual matrix component of a big.matrix is simply a pointer of type <T> (char/short/int/double)

The following code compiles, but the cast to the arma::Mat doesn't work, segfaulting when iterating over the cast matrix.

#include <RcppArmadillo.h>
using namespace Rcpp;

// [[Rcpp::depends(BH, bigmemory, RcppArmadillo)]]
#include <bigmemory/BigMatrix.h>

template <typename T>
void armacast(const arma::Mat<T>& M) {
  // This segfaults
  for (int j = 0; j < 2; j++) {
    for (int i = 0; i < 2; i++) {
      std::cout << M.at(j, i) << std::endl;
    }   
  }
  std::cout << "Success!" << std::endl;                                                                                                                                                                                                       
}

// [[Rcpp::export]]
void armacast(SEXP pDat) {
  XPtr<BigMatrix> xpDat(pDat);

  if (xpDat->matrix_type() == 8) {
    // I can iterate over this *mat and get sensible output.
    double *mat = (double *)xpDat->matrix();
    for (int j = 0; j < 2; j++) {
      for (int i = 0; i < 2; i++) {
        std::cout << *mat + 2 * (j + 0) + i << std::endl;
      }   
    }
    armacast((const arma::Mat<double> &)mat);
  } else {
    std::cout << "Not implemented yet!" << std::endl;
  }
}

In R:

library(Rcpp)
library(RcppArmadillo)
library(bigmemory)
sourceCpp("armacast.cpp")
m <- as.big.matrix(matrix(1:4, 2), type="double")
armacast(m@address)

Answer 1

Great question! We may spin this into another Rcpp Gallery post.

There is one important detail you may have glossed over. Bigmemory objects are external so that we get R to not let its memory management interfere. Armadillo does have constructors for this (and please read the docs and warnings there) so in a first instance we can just do

arma::mat M( (double*) xpDat->matrix(), xpDat->nrow(), xpDat->ncol(), false);

where we using a pointer to the matrix data, as well as row and column counts. A complete version:

// [[Rcpp::export]]
void armacast(SEXP pDat) {
  XPtr<BigMatrix> xpDat(pDat);

  if (xpDat->matrix_type() == 8) {
    arma::mat M(mat, xpDat->nrow(), xpDat>-ncol(), false);
    M.print("Arma matrix M");
  } else {
    std::cout << "Not implemented yet!" << std::endl;
  }
}

It correctly invokes the print method from Armadillo:

R> armacast(m@address)
Arma matrix M
   1.0000   3.0000
   2.0000   4.0000
R>