How to write Rcpp function for simple matrix multiplication in R

Question

I have wrote a Rcpp code to compute element wise matrix multiplication in R. But when try to run this code R stops working and its exiting. How to correct this function?

Thanks in advance.

library(Rcpp)
func <- 'NumericMatrix mmult( NumericMatrix m , NumericMatrix v, bool  byrow=true )
{
if( ! m.nrow() == v.nrow() ) stop("Non-conformable arrays") ;
if( ! m.ncol() == v.ncol() ) stop("Non-conformable arrays") ;

NumericMatrix out(m) ;

for (int i = 1; i <= m.nrow(); i++) 
{
  for (int j = 1; j <= m.ncol(); j++) 
    {
      out(i,j)=m(i,j) * v(i,j) ;
    }
}
return out ;
}'

cppFunction( func )

m1<-matrix(1:4,2,2)
m2<-m1
r1<-mmult(m1,m2)
r2<-m1*m2

With your example, it runs, but `r1` and `r2` are different. — , Jul 28 '15 at 08:28
I want to get my r1 output as r2. How to do it with my code? — kondal, Jul 28 '15 at 08:30
I don't know. Your function works for me with your example, without exiting. Revise your function. — , Jul 28 '15 at 08:32

score 4 · Answer 1 · edited May 23 '17 at 12:32

You have to keep in mind that c++ uses 0 indexed arrays. (See Why does the indexing start with zero in 'C'? and Why are zero-based arrays the norm? .)

So you need to define your loop to run from 0 to m.nrow() - 1

Try this:

func <- '
NumericMatrix mmult( NumericMatrix m , NumericMatrix v, bool  byrow=true )
{
  if( ! m.nrow() == v.nrow() ) stop("Non-conformable arrays") ;
  if( ! m.ncol() == v.ncol() ) stop("Non-conformable arrays") ;

  NumericMatrix out(m) ;

  for (int i = 0; i < m.nrow(); i++) 
  {
  for (int j = 0; j < m.ncol(); j++) 
  {
  out(i,j)=m(i,j) * v(i,j) ;
  }
  }
  return out ;
  }
'

Then I get:

> mmult(m1,m2)
     [,1] [,2]
[1,]    1    9
[2,]    4   16

> m1*m2
     [,1] [,2]
[1,]    1    9
[2,]    4   16

Dirk Eddelbuettel · Accepted Answer · 2015-07-28T14:37:50.723

4

The (at least to me) obvious choice is to use RcppArmadillo:

R> cppFunction("arma::mat matmult(arma::mat A, arma::mat B) { return A % B; }", 
+              depends="RcppArmadillo")
R> m1 <- m2 <- matrix(1:4,2,2)
R> matmult(m1,m2)
     [,1] [,2]
[1,]    1    9
[2,]    4   16
R>

as Armadillo is strongly typed, and has an element-by-element multiplication operator (%) which we use in the one-liner it takes.

edited Jul 28 '15 at 14:37

answered Jul 28 '15 at 11:54

Dirk Eddelbuettel

360,940
56
644
725

How to write Rcpp function for simple matrix multiplication in R

2 Answers2

Linked