Edit: with some clues from Dirk's answer below I worked this out, solution in the body of the question now.
I'm sure this must be documented somewhere but my google skills are failing me.
I'm developing an Rcpp package where I didn't think I would need
dependency on Eigen, so I used NumericVector/Matrix
quite extensively.
However, I now need a Cholesky decomp/solve: I know how to do this
using RcppEigen but I need VectorXd/MatrixXd
's. Say I have a P.S.D.
NumericMatrix, A, and a NumericVector, b. I have tried various variations on the following:
Rcpp::NumericVector b(bb);
Rcpp::NumericMatrix A(AA);
Eigen::Map<Eigen::MatrixXd> A_eigen = as<Eigen::Map<Eigen::MatrixXd> >(A);
Eigen::Map<Eigen::VectorXd> b_eigen = as<Eigen::Map<Eigen::VectorXd> >(b);
Eigen::LLT<Eigen::MatrixXd> llt(A_eigen);
Eigen::VectorXd x = llt.solve(b_eigen);
Rcpp::NumericVector xx(wrap(x));
return xx;
One can then compile this with (denoted here by codeString
) using the inline
package as follows:
f = cxxfunction(signature(AA="matrix",bb="numeric"),codeString,plugin="RcppEigen")
I am aware that I could directly get the Eigen::MatrixXd/VectorXd
from the SEXP instances (e.g. using Eigen::Map) but in my use I need to get these from NumericMatrix/Vector instances.
My A will be pretty small so I'm not too
worried if a full copy gets created. If anyone can offer a clunky
solution, that would be a good, and an elegant solution would be
great!
Thanks in advance, and also for all the great work on Rcpp(Eigen),
David