[Rd] Problem using F77_CALL(dgemm) in a package
Dirk Eddelbuettel
edd at debian.org
Sun Feb 20 16:50:43 CET 2011
On 19 February 2011 at 16:50, Jason Rudy wrote:
 Dear Rdevel,

 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 higherlevel 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 rcppdevel list is open for help and further questions should you have any.
Cheers, Dirk

Dirk Eddelbuettel  edd at debian.org  http://dirk.eddelbuettel.com
More information about the Rdevel
mailing list