I'm writing an R package using Rcpp and RcppEigen and I'm having a problem with matrix slicing and sub-setting. I need to get an off-diagonal square matrix from a larger square matrix
The Eigen::MatrixXd slicing methods seem to be the problem. Using Eigen::seq() method doesn't work at all from R because "no member named 'seq' in namespace 'Eigen'" and the MatrixXd.block(i, j, n, n) method is causing a crash with certain large(ish) values of n. Sometimes it works perfectly fine, but if I increase the size, it causes a fatal crash.
Here's an example of the C++ code:
// [[Rcpp::depends(RcppEigen)]]
#include <iostream>
#include <RcppEigen.h>
using namespace Rcpp;
using Eigen::Map;
using Eigen::MatrixXd;
typedef Map<MatrixXd> MapMatd;
// [[Rcpp::export]]
List crosspart_worker_cpp(const MapMatd& Vij, ...){
int n = Vij.cols()/2;
/* ... some arbitrary code ... */
// extract sub-block of varcovar matrix (only unique pairs)
MatrixXd Vsub = Vij.block(1, n + 1, n, n);
/* ... more code ... */
List out_lst = List::create(Named("Vsub") = Vsub, ...);
return out_lst;
}
and R code that uses it:
# test_crosspart-worker-cpp.R
# setup ----
## source relevant file
# normally build and load package instead of sourceCpp()
if(!exists("crosspart_worker_cpp")){
Rcpp::sourceCpp("crosspart-worker.cpp")
}
## make reproducible
set.seed(75)
## set parameters
n = 100 # rows in original model matrix X
## ... additional setup code ...
# Generate dummy data with arbitrary contents, but correct structure ----
## varcovar matrix
Vij <- matrix(abs(rnorm(2*n * 2*n)), nrow = 2*n)
## ... other code ...
# use the function ----
result <- crosspart_worker_cpp(Vij = Vij, ...)
Why doesn't this work consistently? Are there other options for sub-setting a matrix in RcppEigen?
My original post contained a larger C++ file and I didn't know where the error was. I've been able to identify the problem line of code and the function works perfectly if I remove it. Currently, I am obtaining the matrix subset in R and then passing that as an argument to the C++ function. However, it would be be highly beneficial for my package to do this in all in a C++ function.
I used Rcout to find the problem line of code:
MatrixXd Vsub = VDiag.block(1, n + 1, n, n);
Which is supposed to take the block matrix starting at row 1 and column np and contain np rows and columns according to Eigen: Slicing and Indexing
this is meant to replicate the R code Vsub[1:np, (n + 1):(2*n)].
Why would this cause a crash? Does the Eigen slice method not work in R? is there an RcppEigen specific way? I haven't found one online.
System specs:
OS: Microsoft Windows 10 Education Edition
Processor: AMD Ryzen 7 3700X 8-core
Installed Physical Memory (RAM): 16.0 GB
R version: 4.0.2 (2020-06-22) -- "Taking Off Again"
Rstudio version: 1.3.1056