RcppCore / RcppEigen

Rcpp integration for the Eigen templated linear algebra library
Other
110 stars 40 forks source link

Long vectors are not supported #54

Closed rstub closed 6 years ago

rstub commented 6 years ago

Triggered by this SO question: https://stackoverflow.com/questions/50504929/rcppeigen-cannot-return-matrix-to-r-with-more-than-231-elements

Does it make sense to revisit #10 now that Rcpp supports long vectors? Example:

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

// [[Rcpp::export]]
void get_length_rcpp(Rcpp::IntegerMatrix m){
  Rcpp::Rcout << m.nrow() << ' ' << m.ncol() << ' ' 
              << (m.nrow() * m.ncol()) << ' ' << m.size(); 
}

// [[Rcpp::export]]
void get_length_eigen(Eigen::Map<Eigen::MatrixXi> m){
  Rcpp::Rcout << m.rows() << ' ' << m.cols() << ' ' 
              << (m.rows() * m.cols()) << ' ' << m.size(); 
}

/*** R
N <- 5e4
A <- matrix(1L, ncol = N, nrow = N)
get_length_rcpp(A)
get_length_eigen(A)
*/

Output:

> N <- 50000

> A <- matrix(1, ncol = N, nrow = N)

> get_length_rcpp(A)
50000 50000 -1794967296 2500000000
> get_length_eigen(A)
Error in get_length_eigen(A) : 
  long vectors not supported yet: ../../src/include/Rinlinedfuns.h:519
Calls: <Anonymous> ... withVisible -> eval -> eval -> get_length_eigen -> .Call
Execution halted

I am using RcppEigen 0.3.3.4.0 with R 3.5.0.

eddelbuettel commented 6 years ago

A careful pull request would probably be looked upon favorably. You can look at @thirdwing 's work for these changes in Rcpp itself; this should be doable along the same lines.

rstub commented 6 years ago

Ok, I am on it. Here is a second test case that would not be covered by #10:

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

// [[Rcpp::export]]
SEXP create_matrix(R_xlen_t n) {
    Eigen::MatrixXi A = Eigen::MatrixXi::Zero(n, n);
    return Rcpp::wrap(A);
}

/*** R
N <- 5e4
length(create_matrix(N))
*/

Error message:

Error in create_matrix(N) : negative length vectors are not allowed

eddelbuettel commented 6 years ago

But isn't R_xlen_t really a double passing it to the Eigen ctor there is probably not a good idea. What if you cast that to, say, a int64_t or whatever Eigen uses?

(Edit: Nope, I misremembered. bigmemory used the double trick, R is indeed ptrdiff_t as you show below.)

rstub commented 6 years ago

It is of type ptrdiff_t:

#ifdef LONG_VECTOR_SUPPORT
    typedef ptrdiff_t R_xlen_t;
# define R_XLEN_T_MAX 4503599627370496
# define R_SHORT_LEN_MAX 2147483647
#else
    typedef int R_xlen_t;
# define R_XLEN_T_MAX R_LEN_T_MAX
#endif

R_XLEN_T_MAX equals 2^52. Probably that is used to check that long vectors can be indexed via doubles.