[Rcpp-devel] Fwd: Element-wise matrix operations

[Jason Nielsen]

Hi all,

I was wondering if there is a technical reason that element-wise matrix
operations between Rcpp's Matrix type aren't implemented. There is
vector/vector, vector/scalar and matrix/scalar but no matrix/matrix. The
following naive code:

#include <Rcpp.h>
using namespace Rcpp;

NumericMatrix operator+(NumericMatrix A, NumericMatrix B) {
if ((A.nrow() != B.nrow()) | (A.ncol() != B.ncol()))
stop("Dimensions of the two matrices do not conform!");
NumericVector vA = as<NumericVector>(A);
NumericVector vB = as<NumericVector>(B);
NumericVector vC = vA + vB;
NumericMatrix C = NumericMatrix(A.nrow(), A.ncol(), vC.begin());
return C;
}

// [[Rcpp::export]]
NumericMatrix MatAdd(const NumericMatrix A, const NumericMatrix B) {
return A+B;
}

works. From my brief fiddling it seems there is no copying going on above
so I don't think there is too much extra overhead in the abstraction. My
guess is that it just hasn't been a priority and "patches happily accepted"
is the answer but I'm curious if there is some other less obvious reason as
it seems to me all the heavy lifting was done to implement vector/vector.

Many thanks,
Jason

[Dirk Eddelbuettel]

On 20 January 2018 at 17:45, Jason Nielsen wrote:
| Hi all,
|
| I was wondering if there is a technical reason that element-wise matrix
| operations between Rcpp's Matrix type aren't implemented. There is
| vector/vector, vector/scalar and matrix/scalar but no matrix/matrix. The
| following naive code:
|
| #include <Rcpp.h>
| using namespace Rcpp;
|
| NumericMatrix operator+(NumericMatrix A, NumericMatrix B) {
| if ((A.nrow() != B.nrow()) | (A.ncol() != B.ncol()))
| stop("Dimensions of the two matrices do not conform!");
| NumericVector vA = as<NumericVector>(A);
| NumericVector vB = as<NumericVector>(B);
| NumericVector vC = vA + vB;
| NumericMatrix C = NumericMatrix(A.nrow(), A.ncol(), vC.begin());
| return C;
| }
|
| // [[Rcpp::export]]
| NumericMatrix MatAdd(const NumericMatrix A, const NumericMatrix B) {
| return A+B;
| }
|
| works. From my brief fiddling it seems there is no copying going on above
| so I don't think there is too much extra overhead in the abstraction. My
| guess is that it just hasn't been a priority and "patches happily accepted"
| is the answer but I'm curious if there is some other less obvious reason as
| it seems to me all the heavy lifting was done to implement vector/vector.

R code:

sample(c("History", "Lazyness"), 1) ## just kidding. or maybe not.

Years back I usually pointed at RcppArmadillo if you wanted "math with
matrices" as we just didn't have all the operators. Now, we have been filling
them in, so this operator+() makes some sense. Doing them "well" in the
template code is a bit more involved but mayne you want to poke around?

I am sure some of the regulars here would be happy to help you along.

Dirk

[Jason Nielsen]

Thank you for your quick response. I was only thinking about the
element-wise matrix/matrix equivalents that already have vector/vector
implementations. My first look was telling me that this shouldn't be
difficult. That said reading and using templated C++ code is a different
story than modifying it!

Jason

Post by Dirk Eddelbuettel
| Hi all,
|
| I was wondering if there is a technical reason that element-wise matrix
| operations between Rcpp's Matrix type aren't implemented. There is
| vector/vector, vector/scalar and matrix/scalar but no matrix/matrix. The
|
| #include <Rcpp.h>
| using namespace Rcpp;
|
| NumericMatrix operator+(NumericMatrix A, NumericMatrix B) {
| if ((A.nrow() != B.nrow()) | (A.ncol() != B.ncol()))
| stop("Dimensions of the two matrices do not conform!");
| NumericVector vA = as<NumericVector>(A);
| NumericVector vB = as<NumericVector>(B);
| NumericVector vC = vA + vB;
| NumericMatrix C = NumericMatrix(A.nrow(), A.ncol(), vC.begin());
| return C;
| }
|
| // [[Rcpp::export]]
| NumericMatrix MatAdd(const NumericMatrix A, const NumericMatrix B) {
| return A+B;
| }
|
| works. From my brief fiddling it seems there is no copying going on above
| so I don't think there is too much extra overhead in the abstraction. My
| guess is that it just hasn't been a priority and "patches happily accepted"
| is the answer but I'm curious if there is some other less obvious reason as
| it seems to me all the heavy lifting was done to implement vector/vector.
sample(c("History", "Lazyness"), 1) ## just kidding. or maybe not.
Years back I usually pointed at RcppArmadillo if you wanted "math with
matrices" as we just didn't have all the operators. Now, we have been filling
them in, so this operator+() makes some sense. Doing them "well" in the
template code is a bit more involved but mayne you want to poke around?
I am sure some of the regulars here would be happy to help you along.
Dirk
--