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?

Improve this question
6
  • 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"). Commented Dec 1, 2021 at 1:40
  • 1
    I allowed myself to edit your post title: this had nothing to do with 'execute'. Commented Dec 1, 2021 at 1:41
  • thanks a lot @DirkEddelbuettel for that suggestion and for responding! Let me give it a try Commented Dec 1, 2021 at 1: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 Commented Dec 1, 2021 at 2: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 ... Commented Dec 1, 2021 at 4:26

1 Answer 1

Reset to default
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
Improve this answer
Sign up to request clarification or add additional context in comments.

Comments

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.