Mihai
Mihai

Reputation: 2967

Update Rcpp::NumericMatrix passed by reference using RcppArmadillo submat()

Following on this question, I am trying to understand how to efficiently update a subset of a Rccp::NumericMatrix data type.

I have the following scenario:

The code looks like this:

#include <iostream>
#include <RcppArmadillo.h>
// [[Rcpp::depends(RcppArmadillo)]]


// [[Rcpp::export]]
void updateMatrix(const Rcpp::NumericMatrix &m)
{
    std::cout << m << std::endl;

    Rcpp::as<arma::mat>(m).submat(0, 0, 3, 3) = Rcpp::as<arma::mat>(m).submat(0, 0, 3, 3) + 1;

    std::cout << m << std::endl;
}

To run it from R I use:

m = matrix(0, 5, 5)

updateMatrix(m)

And the results are:

> updateMatrix(m)
0.00000 0.00000 0.00000 0.00000 0.00000
0.00000 0.00000 0.00000 0.00000 0.00000
0.00000 0.00000 0.00000 0.00000 0.00000
0.00000 0.00000 0.00000 0.00000 0.00000
0.00000 0.00000 0.00000 0.00000 0.00000

0.00000 0.00000 0.00000 0.00000 0.00000
0.00000 0.00000 0.00000 0.00000 0.00000
0.00000 0.00000 0.00000 0.00000 0.00000
0.00000 0.00000 0.00000 0.00000 0.00000
0.00000 0.00000 0.00000 0.00000 0.00000

It's the first time I'm using Rcpp and RcppArmadillo and they are absolutely amazing. I appreciate any help with this scenario.

Upvotes: 0

Views: 1096

Answers (1)

Bob Jansen
Bob Jansen

Reputation: 1287

The left side of your assignment in updateMatrix creates a temporary that is discarded after assignment. Therefore, m doesn't change at all. The code can't work as you expected as the that would mean the type of m would change. Look below:

#include <typeinfo>
#include <iostream>
#include <RcppArmadillo.h>
// [[Rcpp::depends(RcppArmadillo)]]


// [[Rcpp::export]]
void updateMatrix(const Rcpp::NumericMatrix &m)
{
  std::cout << m << std::endl;

  std::cout << typeid(m).name() << std::endl;

  arma::mat m2 = Rcpp::as<arma::mat>(m);

  std::cout << typeid(m2).name() << std::endl;

  m2.submat(0, 0, 3, 3) = Rcpp::as<arma::mat>(m).submat(0, 0, 3, 3) + 1;

  std::cout << m2 << std::endl;
}

Running this code gives:

> m = matrix(0, 5, 5)
> updateMatrix(m)
0.00000 0.00000 0.00000 0.00000 0.00000
0.00000 0.00000 0.00000 0.00000 0.00000
0.00000 0.00000 0.00000 0.00000 0.00000
0.00000 0.00000 0.00000 0.00000 0.00000
0.00000 0.00000 0.00000 0.00000 0.00000

N4Rcpp6MatrixILi14ENS_15PreserveStorageEEE
N4arma3MatIdEE
   1.0000   1.0000   1.0000   1.0000        0
   1.0000   1.0000   1.0000   1.0000        0
   1.0000   1.0000   1.0000   1.0000        0
   1.0000   1.0000   1.0000   1.0000        0
        0        0        0        0        0

Upvotes: 3

Related Questions