On 19 February 2011 at 16:50, Jason Rudy wrote:
| Dear R-devel,
|
| I've written a numerical solver for SOCPs (second order cone programs)
| in R, and now I want to move most of the solver code into C for speed.
| I've written combined R/C packages before, but in this case I need to
| do matrix operations in my C code. As I have never done that before,
| I'm trying to write some simple examples to make sure I understand the
| basics. I am stuck on the first one. I'm trying to write a function
| to multiply two matrices using the blas routine dgemm. The name of my
| example package is CMATRIX. My code is as follows.
There is of course merit in working through the barebones API but in case you
would consider a higher-level alternative, consider these few lines based on
RcppArmadillo (which end up calling dgemm() for you via R's linkage to the BLAS)
suppressMessages(library(inline))
txt <- '
Rcpp::NumericMatrix Ar(A); // creates Rcpp matrix from SEXP
Rcpp::NumericMatrix Br(B); // creates Rcpp matrix from SEXP
arma::mat Am(Ar.begin(), Ar.nrow(), Ar.ncol(), false); // Arma mat, no copy
arma::mat Bm(Br.begin(), Br.nrow(), Br.ncol(), false); // Arma mat, no copy
arma::mat Cm = Am * Bm;
return Rcpp::wrap(Cm);
'
mmult <- cxxfunction(signature(A="numeric", B="numeric"),
body=txt,
plugin="RcppArmadillo")
A <- matrix(1:9, 3, 3)
B <- matrix(9:1, 3, 3)
C <- mmult(A, B)
print(C)
which when passed into R yield
R> txt <- '
+ Rcpp::NumericMatrix Ar(A); // creates Rcpp matrix from SEXP
+ Rcpp::NumericMatrix Br(B); // creates Rcpp matrix from SEXP
+ arma::mat Am(Ar.begin(), Ar.nrow(), Ar.ncol(), false); // Arma mat, no
copy
+ arma::mat Bm(Br.begin(), Br.nrow(), Br.ncol(), false); // Arma mat, no
copy
+ arma::mat Cm = Am * Bm;
+ return Rcpp::wrap(Cm);
+ '
R> mmult <- cxxfunction(signature(A="numeric", B="numeric"),
+ body=txt,
+ plugin="RcppArmadillo")
R> A <- matrix(1:9, 3, 3)
R> B <- matrix(9:1, 3, 3)
R> C <- mmult(A, B)
R> C
[,1] [,2] [,3]
[1,] 90 54 18
[2,] 114 69 24
[3,] 138 84 30
R>
You can then use this helper function to have a package created for you:
RcppArmadillo.package.skeleton("mmult", mmult, path="/tmp")
The rcpp-devel list is open for help and further questions should you have any.
Cheers, Dirk
--
Dirk Eddelbuettel | e...@debian.org | http://dirk.eddelbuettel.com
______________________________________________
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
