[Rd] Problem using F77_CALL(dgemm) in a package

[Dirk Eddelbuettel]

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 | edd at debian.org | http://dirk.eddelbuettel.com