1

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?

coatless
  • 20,011
  • 13
  • 69
  • 84
user1993
  • 498
  • 1
  • 10
  • 22
  • 3
    Maybe try simpler Eigen code? You have an complex decomposition there that you expect to be transfered in passing to a list element, then a list and then (finally) a `SEXP` that R can take. Sometimes ... you need to decompose that. As the saying goes: "try to walk before you run". Which is write I write in the Rcpp documentation to first try `Rcpp::evalCpp("2 + 2")`. – Dirk Eddelbuettel Dec 01 '21 at 01:40
  • 1
    I allowed myself to edit your post title: this had nothing to do with 'execute'. – Dirk Eddelbuettel Dec 01 '21 at 01:41
  • thanks a lot @DirkEddelbuettel for that suggestion and for responding! Let me give it a try – user1993 Dec 01 '21 at 01:54
  • 1
    @DirkEddelbuettel somehow that did the job! Thanks a lot. I have provided an answer to this question for anyone who might struggle with the same – user1993 Dec 01 '21 at 02:58
  • 1
    There is no 'somehow'. That is just how it works: when you force six steps into one, it may fail. So decompose one by one ... – Dirk Eddelbuettel Dec 01 '21 at 04:26
  • Gotcha! thanks a lot – user1993 Dec 01 '21 at 04:40

1 Answers1

1

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
user1993
  • 498
  • 1
  • 10
  • 22