Rcpp installed but compilation error from complex code snippet

105 Views Asked by At

I have Rcpp and RccpEigen already installed in RStudio. I am able to run an Rcpp code (that didn't use RccpEigen) successfully as well. However the following code which uses both doesn't seem to work.

Here is the code -

library(Rcpp)
library(RcppEigen)
sourceCpp(code = '
  #include <Rcpp.h>
  #include <RcppEigen.h>
  // [[Rcpp::depends(RcppEigen)]]
  using namespace Rcpp;
  using namespace Eigen;
  using namespace RcppEigen;
  // [[Rcpp::export]]
    List luEigen(MatrixXd M) {
    FullPivLU<MatrixXd> luE(M);
    return List::create(Named("L_matrix") = luE.matrixLU().triangularView<Upper>());
}')

A <- 0.8 + 0.2 * diag(100)
(luEigen(A))

This code gives a really long error, so here are the key error lines -

/Library/Frameworks/R.framework/Versions/4.1/Resources/library/Rcpp/include/Rcpp/generated/Vector__create.h:71:10: note: in instantiation of function template specialization 'Rcpp::Vector<19, PreserveStorage>::create__dispatch<Rcpp::traits::named_object<Eigen::TriangularView<const Eigen::Matrix<double, -1, -1, 0>, 2>>>' requested here
                return create__dispatch( typename traits::integral_constant<bool,
                       ^
file16bbd8305f5c.cpp:11:18: note: in instantiation of function template specialization 'Rcpp::Vector<19, PreserveStorage>::create<Rcpp::traits::named_object<Eigen::TriangularView<const Eigen::Matrix<double, -1, -1, 0>, 2>>>' requested here
    return List::create(Named("L_matrix") = luE.matrixLU().triangularView<Upper>());
                 ^
18 warnings and 1 error generated.
make: *** [file16bbd8305f5c.o] Error 1
clang++ -mmacosx-version-min=10.13 -std=gnu++14 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG   -I"/Library/Frameworks/R.framework/Versions/4.1/Resources/library/Rcpp/include" -I"/Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppEigen/include" -I"/private/var/folders/_3/wdql3v5d4vggzffw3xdcr3p80000gn/T/RtmpQioi38/sourceCpp-x86_64-apple-darwin17.0-1.0.7" -I/usr/local/include   -fPIC  -Wall -g -O2  -c file16bbd8305f5c.cpp -o file16bbd8305f5c.o

Given that Rcpp and RccpEigen are installed and a different Rccp code does work, what may be causing error in this code?

1

There are 1 best solutions below

0
On

With a very helpful suggestion from @Dirk, I simplified the decomposition and that did the trick. Still not sure why the more complex construction threw an error, but bottom line is that simplification gets the job done. Here is my modified code that works -

library(Rcpp)
library(RcppEigen)
sourceCpp(code = '
  #include <Rcpp.h>
  #include <RcppEigen.h>
  // [[Rcpp::depends(RcppEigen)]]
  using namespace Rcpp;
  using namespace Eigen;
  using namespace RcppEigen;
  // [[Rcpp::export]]
    List luEigen(MatrixXd M) { // here I name our function 
    FullPivLU<MatrixXd> luE(M); // here I perform the decomposition
    MatrixXd upper = luE.matrixLU().triangularView<Upper>(); // this creates the upper matrix
    MatrixXd lower = luE.matrixLU().triangularView<StrictlyLower>(); // this creates the lower matrix
    return List::create(Named("U_matrix") = upper, Named("L_matrix") = lower); // this makes the list of the 2 matrices
}')

A <- 0.8 + 0.2 * diag(100)
(luEigen(A))

You could possibly speed it up further by doing the decomposition only once and calling the upper and lower triangular matrices from it, like so -

FullPivLU<MatrixXd> luE(M); // here I perform the decomposition
MatrixXd decomp = luE.matrixLU();
MatrixXd upper = decomp.triangularView<Upper>(); // this creates the upper matrix
MatrixXd lower = decomp.triangularView<StrictlyLower>(); // this creates the lower matrix