如何在 R 中编写用于简单矩阵乘法的 Rcpp 函数

Question

我已经编写了一个 Rcpp 代码来计算 R 中的逐元素矩阵乘法。但是当尝试运行时，此代码 R 停止工作并退出。如何修正这个功能？

提前致谢。

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

Answer 1

你必须记住c++ uses 0 indexed arrays. (See Why does the indexing start with zero in 'C'? and Why are zero-based arrays the norm?。）

因此您需要将循环定义为运行从 0 到 m.nrow() - 1

试试这个：

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 ;
  }
'

然后我得到：

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

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

Answer 2

（至少对我而言）显而易见的 选择是使用 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>

因为犰狳是强类型的，并且有一个逐个元素的乘法运算符（%），我们在它需要的一行中使用它。

如何在 R 中编写用于简单矩阵乘法的 Rcpp 函数

How to write Rcpp function for simple matrix multiplication in R

r

rcpp