5

I'm working on a package that requires some very fast matrix multiplication so looking to use RcppEigen. For a variety of reasons though having to do with the need for multidimensional arrays, I need to convert a created object of class Eigen::MatrixXd to class Rcpp::NumericMatrix.

I tried reversing the steps listed in RcppEigen::FastLm.cpp, but that doesn't seem to work

e.g. instead of using

const Map<MatrixXd> X(as<Map<MatrixXd> >(Xs));

I tried

Rcpp:NumericMatrix X(as<Rcpp::NumericMatrix>(Xs));

where Xs is a matrix of class Eigen::MatrixXd but that didn't seem to work:" error: no matching function for call to 'as' return Rcpp::asRcpp::NumericMatrix(z);"

If this isn't at all possible I can try another direction.

Basically what I need to do in R speak is

a = matrix(1, nrow = 10, ncol = 10)

b = array(0, c(10,10,10))

b[,,1] = a

To give a clearer starting example

How would I go about storing an object of class MatrixXd in an object of class NumericMatrix?

#include <Rcpp.h>
#include <RcppEigen.h>
using namespace Rcpp;
using namespace Eigen;

// [[Rcpp::export]]
NumericMatrix sample_problem() {

  Eigen::MatrixXd x(2, 2); x << 1,1,2,2;

  Eigen::MatrixXd z(2, 2);

  Eigen::MatrixXd y(2,2); y << 3,3,4,4;

  z =  x * y; // do some eigen matrix multiplication

  Rcpp::NumericMatrix w(2,2);

  // what I'd like to be able to do somehow: 
  // store the results of the eigen object z in
  // a NumericMatrix w
  // w = z;

  return w;
} 

DanO
  • 600
  • 3
  • 11
  • 1
    Can you focus and sharpen the question and actually show code? Ex ante a _multidimensional array_ simply won't fit into a 2d NumericMatrix so maybe you can make it more concrete about 2-d slices of a multi-d array? And maybe once you have 2-d matrix slices or view in Eigen it will become more obvious how those become a SEXP which can of course become a NumericMatrix (and all can happen zero-copy). – Dirk Eddelbuettel Jun 26 '20 at 03:45
  • Thanks, see hopefully sharpened question with code added to original post – DanO Jun 26 '20 at 15:03
  • Perfect, see below. – Dirk Eddelbuettel Jun 26 '20 at 15:14

1 Answers1

5

Thanks for posting code! It makes everything easier. I just rearranged you code the tiniest bit.

The key changes is to "explicitly" go back from the Eigen representation via an RcppEigen helper to a SEXP, and to then instantiate the matrix. Sometimes ... the compiler needs a little nudge.

Code

#include <Rcpp.h>
#include <RcppEigen.h>

// [[Rcpp::depends(RcppEigen)]]

// [[Rcpp::export]]
Rcpp::NumericMatrix sample_problem() {

  Eigen::MatrixXd x(2, 2), y(2, 2);
  x << 1,1,2,2;
  y << 3,3,4,4;

  // do some eigen matrix multiplication
  Eigen::MatrixXd z =  x * y;

  // what I'd like to be able to do somehow:
  // store the results of the eigen object z in
  // a NumericMatrix w
  // w = z;
  SEXP s = Rcpp::wrap(z);
  Rcpp::NumericMatrix w(s);

  return w;
}

/*** R
sample_problem()
*/

Demo

R> sourceCpp("demo.cpp)

R> sample_problem()
     [,1] [,2]
[1,]    7    7
[2,]   14   14
R> 

With g++-9 I need -Wno-ignored-attributes or I get screens and screens of warnings...

Dirk Eddelbuettel
  • 360,940
  • 56
  • 644
  • 725
  • That's exactly what I was looking for, thanks! Now to tackle assignment to a multidimensional array, but this should get me going in the right direction – DanO Jun 26 '20 at 15:37
  • Slight issue with this I think: running into an error "Not compatible with requested type: [type=any; target=double]." seemingly at random; running the exact same function again after encountering that error works just fine. There are no stochastic processes in the model. Seems potentially related to this issue? https://stackoverflow.com/questions/32375537/error-not-compatible-with-requested-type-happening-at-random. Apologies that this is a vague question but I am stumped; any general thoughts/tips are appreciated! – DanO May 07 '21 at 18:38
  • No idea, doing `Rcpp::sourceCpp("~/git/stackoverflow/62586950/answer.cpp")` still works fine here (modulo the lots of Eigen compilation noise). – Dirk Eddelbuettel May 07 '21 at 19:00
  • Sorry, should have clarified, the error is in the code I built based on this suggestion (https://github.com/DanOvando/marlin/blob/e19311d387a9dedeea20794df27f1e43ab5fc465/src/sim-fish.cpp), but guessing if it's not an obvious and known problem to you it's something very weird I need to do more digging on. thanks – DanO May 07 '21 at 19:10