RcppEigen slicing method to obtain subset of a matrix not working consistently and causes fatal error/crash

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

Recommended topics

Hot tags