I'm new to Rcpp and Eigen and I'm trying to run a inplace matrix multiplication while I run into this error. R didn't tell me what the error is about.
SEXP cpp_hom_crit(
const Eigen::Map<Eigen::MatrixXd> Q, Eigen::Map<Eigen::VectorXd> d, Eigen::Map<Eigen::MatrixXd> Qinv) {
Q *= d.asDiagonal() * Qinv; //error happens here
return Rcpp::wrap(Q);
}
I tried this one also, but it does not work either
SEXP cpp_hom_crit(
const Eigen::Map<Eigen::MatrixXd> Q, Eigen::Map<Eigen::VectorXd> d, Eigen::Map<Eigen::MatrixXd> Qinv) {
Q = Q * d.asDiagonal() * Qinv; //error happens here
return Rcpp::wrap(Q);
}
How can I resolve this error and possibly make this faster?
Here is one way to fix your problem. As @chtz alluded to in the comment, you cannot declare a variable const and assign to it. So here we assign to a new temporary variable and return that result.
#include <RcppEigen.h>
// [[Rcpp::depends(RcppEigen)]]
// [[Rcpp::export]]
Eigen::MatrixXd cpp_hom_crit(const Eigen::Map<Eigen::MatrixXd> Q,
const Eigen::Map<Eigen::VectorXd> d,
const Eigen::Map<Eigen::MatrixXd> Qinv) {
auto res = Q * d.asDiagonal() * Qinv;
return res;
}
/*** R
v <- as.numeric(1:3)
M <- v %*% t(v)
Mi <- chol2inv(M)
cpp_hom_crit(M, v, Mi)
*/
> Rcpp::sourceCpp("~/git/stackoverflow/76625654/answer.cpp")
> v <- as.numeric(1:3)
> M <- v %*% t(v)
> Mi <- chol2inv(M)
> cpp_hom_crit(M, v, Mi)
[,1] [,2] [,3]
[1,] 0.75 0.0694444 0.0370370
[2,] 1.50 0.1388889 0.0740741
[3,] 2.25 0.2083333 0.1111111
>
The values are of course nonsensical but they illustrate at least the correct dimension. Note that libraries like LAPACK often have dedicated functions for common vector-matrix product expressions, do the looping internally and may also avoid the temporary so not it is not clear that Eigen will automatically be (much) faster.