cmo
cmo

Reputation: 4114

RcppEigen matrix multiplcation only returns first element, repeated

Simple matrix multiplication in RcppEigen gives a matrix of correct dimensions, but the first element is repeated to all the elements.

Rcpp source code using RcppEigen for matrix multiply:

#include <Rcpp.h>
#include <RcppEigen.h>


using namespace Rcpp;

// [[Rcpp::depends(RcppEigen)]]
// [[Rcpp::export]]
NumericMatrix myfun(const NumericVector a ,  const NumericVector b) {


        const Eigen::Map<Eigen::MatrixXd> a_eig(Rcpp::as<Eigen::Map<Eigen::MatrixXd> >(a));
        const Eigen::Map<Eigen::MatrixXd> b_eig(Rcpp::as<Eigen::Map<Eigen::MatrixXd> >(b));

        return wrap(a_eig * b_eig);
}

Calling from R:

a=matrix(data=1,nrow=10,ncol=1)
b=sample(100,10)
a * b


      [,1]
 [1,]   67
 [2,]   59
 [3,]   19
 [4,]   68
 [5,]   83
 [6,]    4
 [7,]   28
 [8,]   88
 [9,]   97
[10,]   43

myfun(a,b)

      [,1]
 [1,]   67
 [2,]   67
 [3,]   67
 [4,]   67
 [5,]   67
 [6,]   67
 [7,]   67
 [8,]   67
 [9,]   67
[10,]   67

Upvotes: 2

Views: 138

Answers (1)

Dirk is no longer here
Dirk is no longer here

Reputation: 368509

You are going too fast. Use intermediate results -- your product is now what you think it is. And even with the rest of what you have, you are getting confused about when Eigen::Map makes sense and when not. You also do not need as<> and wrap() -- RcppEigen takes care of that.

Code

#include <Rcpp.h>
#include <RcppEigen.h>

// [[Rcpp::depends(RcppEigen)]]

// [[Rcpp::export]]
Eigen::MatrixXd myfun(Eigen::Map<Eigen::MatrixXd> a,
                      Eigen::Map<Eigen::MatrixXd> b) {
  Eigen::MatrixXd res = a * b;
  return res;
}

Output

 R> Rcpp::sourceCpp("~/git/stackoverflow/55797031/answer.cpp")
 R> myfun(matrix(2,2,2), matrix(3,2,2))
      [,1] [,2]
 [1,]   12   12
 [2,]   12   12
 R>
 R> myfun(matrix(as.numeric(1:4),2,2), matrix(as.numeric(4:1),2,2))
      [,1] [,2]
 [1,]   13    5
 [2,]   20    8
 R>

Upvotes: 4

Related Questions