2015-07-28 2 views
1

Я написал код Rcpp для вычисления умножения элемента в матрице в R. Но при попытке запустить этот код R перестает работать и выходит. Как исправить эту функцию?Как написать функцию Rcpp для простого умножения матрицы в 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 
+0

В вашем примере это работает, но 'r1' и' r2' отличаются. –

+0

Я хочу получить свой выход r1 как r2. Как это сделать с моим кодом? – kondal

+0

Не знаю. Ваша функция работает для меня с вашим примером, без выхода. Пересмотрите свою функцию. –

ответ

3

В (по крайней мере, для меня) очевидно выбор должен использовать 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> 

, как Armadillo сильно типизированный и имеет элемент-поэлементный оператор умножения (%), который мы используем в однострочнике он принимает.

4

Вы должны иметь в виду, что использует 0 индексированных массивов. (См Why does the indexing start with zero in 'C'? и 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