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.
I have a file matrix.c in my src directory:
#include <R.h>
#include <R_ext/Utils.h>
#include <R_ext/Lapack.h>
#include <R_ext/BLAS.h>
//Computes C = A*B
void R_matrix_multiply(double * A, double * B, int * m, int *n, int *
p, double * C){
double one = 1.0;
double zero = 0.0;
//Just printing the input arguments
Rprintf("m = %d, n = %d, p = %d\n",*m,*n,*p);
int i;
for(i=0;i<(*m**n);i++){
Rprintf("%f ",A[i]);
}
Rprintf("\n");
for(i=0;i<(*n**p);i++){
Rprintf("%f ",B[i]);
}
Rprintf("\n");
for(i=0;i<(*m**p);i++){
Rprintf("%f ",C[i]);
}
Rprintf("\n");
//Here is the actual multiplication
F77_CALL(dgemm)("N","N",m,n,p,&one,A,m,B,n,&zero,C,m);
}
And the file C_matrix_multiply.R in my R directory:
C_matrix_multiply = function(A,B){
C <- matrix(0,nrow(A),ncol(B))
cout <-
.C("R_matrix_multiply",as.double(A),as.double(B),nrow(A),ncol(A),ncol(B),as.double(C))
return(matrix(cout$C,nrowA,ncol(B)))
}
My namespace file is:
export("C_matrix_multiply")
useDynLib(CMATRIX.so,R_matrix_multiply)
I'm not sure if it's necessary, but I've also included a Makevars.in
file in my src directory:
PKG_CPPFLAGS=@PKG_CPPFLAGS@
PKG_CFLAGS=@PKG_CFLAGS@
PKG_LIBS=@PKG_LIBS@ ${LAPACK_LIBS} ${BLAS_LIBS} ${FLIBS}
which I simply copied from the diversitree package, which seems to use
a lot of fortran. I have the same problem (which I am getting to)
with or without this Makevars.in file.
I install my package using:
R CMD INSTALL CMATRIX
Then I start up R and attempt to run the following code:
#Make some random matrices
A = matrix(rnorm(8),4,2)
B = matrix(rnorm(6),2,3)
#Load my package
library(CMATRIX)
#Print the matrices
A
B
#Try to multiply them
product = C_matrix_multiply(A,B)
What I want, and what according to my understanding should happen, is
for product to contain the same matrix as would result from A %*% B.
Instead, I get the following:
> A = matrix(rnorm(8),4,2)
> B = matrix(rnorm(6),2,3)
> library(CMATRIX)
> A
[,1] [,2]
[1,] -0.4981664 -0.7243532
[2,] 0.1428766 -1.5501623
[3,] -2.0624701 1.5104507
[4,] -0.5871962 0.3049442> B
[,1] [,2] [,3]
[1,] 0.02477964 0.5827084 1.8434375
[2,] -0.20200104 1.7294264 0.9071397> C_matrix_multiply(A,B)
m = 4, n = 2, p = 3
-0.498166 0.142877 -2.062470 -0.587196 -0.724353 -1.550162 1.510451 0.304944
0.024780 -0.202001 0.582708 1.729426 1.843437 0.907140
0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000
0.000000 0.000000 0.000000 0.000000 0.000000
Parameter 10 to routine DGEMM was incorrect
Mac OS BLAS parameter error in DGEMM , parameter #0, (unavailable), is 0
and R immediately dies. I know the arguments are being passed into
the C code and everything up to my F77_CALL is functioning based on
the printed output. The problem is definitely something to do with my
F77_CALL(dgemm) line. My understanding is that parameter 10 should be
the leading dimension of the matrix B, which in this case should be
equal to 2, the number of rows in that matrix, which is what I am
doing. I have also considered that parameter numbering starts at 0,
in which case the incorrect parameter is &zero, but again that seems
correct to me. All of my reading and research suggests I am doing
everything correctly, so I am somewhat stumped. Perhaps I am missing
something simple or obvious, as I have never done this before and am
proceeding with only google and the R docs as my guide. I am
wondering if anybody can see what I'm doing wrong here, or perhaps
something I could do to try to fix it. Any assistance would be
greatly appreciated.
Best Regards,
Jason Rudy
Graduate Student
Bioinformatics and Medical Informatics Program
San Diego State University
Look a close look at matprod in src/main/array in the R sources. Hint: it is the other dimensions you have wrong. And as BLAS is Fortran, counts do start at 1. On Sat, 19 Feb 2011, 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. > > I have a file matrix.c in my src directory: > > #include <R.h> > #include <R_ext/Utils.h> > #include <R_ext/Lapack.h> > #include <R_ext/BLAS.h> > > //Computes C = A*B > void R_matrix_multiply(double * A, double * B, int * m, int *n, int * > p, double * C){ > double one = 1.0; > double zero = 0.0; > > //Just printing the input arguments > Rprintf("m = %d, n = %d, p = %d\n",*m,*n,*p); > int i; > for(i=0;i<(*m**n);i++){ > Rprintf("%f ",A[i]); > } > Rprintf("\n"); > for(i=0;i<(*n**p);i++){ > Rprintf("%f ",B[i]); > } > Rprintf("\n"); > for(i=0;i<(*m**p);i++){ > Rprintf("%f ",C[i]); > } > Rprintf("\n"); > > > //Here is the actual multiplication > F77_CALL(dgemm)("N","N",m,n,p,&one,A,m,B,n,&zero,C,m); > } > > And the file C_matrix_multiply.R in my R directory: > > C_matrix_multiply = function(A,B){ > C <- matrix(0,nrow(A),ncol(B)) > cout <- .C("R_matrix_multiply",as.double(A),as.double(B),nrow(A),ncol(A),ncol(B),as.double(C)) > return(matrix(cout$C,nrowA,ncol(B))) > > } > > My namespace file is: > > export("C_matrix_multiply") > useDynLib(CMATRIX.so,R_matrix_multiply) > > I'm not sure if it's necessary, but I've also included a Makevars.in > file in my src directory: > > PKG_CPPFLAGS=@PKG_CPPFLAGS@ > PKG_CFLAGS=@PKG_CFLAGS@ > PKG_LIBS=@PKG_LIBS@ ${LAPACK_LIBS} ${BLAS_LIBS} ${FLIBS} > > which I simply copied from the diversitree package, which seems to use > a lot of fortran. I have the same problem (which I am getting to) > with or without this Makevars.in file. > > I install my package using: > > R CMD INSTALL CMATRIX > > Then I start up R and attempt to run the following code: > > #Make some random matrices > A = matrix(rnorm(8),4,2) > B = matrix(rnorm(6),2,3) > > #Load my package > library(CMATRIX) > > #Print the matrices > A > B > > #Try to multiply them > product = C_matrix_multiply(A,B) > > What I want, and what according to my understanding should happen, is > for product to contain the same matrix as would result from A %*% B. > Instead, I get the following: > >> A = matrix(rnorm(8),4,2) >> B = matrix(rnorm(6),2,3) >> library(CMATRIX) >> A > [,1] [,2] > [1,] -0.4981664 -0.7243532 > [2,] 0.1428766 -1.5501623 > [3,] -2.0624701 1.5104507 > [4,] -0.5871962 0.3049442 >> B > [,1] [,2] [,3] > [1,] 0.02477964 0.5827084 1.8434375 > [2,] -0.20200104 1.7294264 0.9071397 >> C_matrix_multiply(A,B) > m = 4, n = 2, p = 3 > -0.498166 0.142877 -2.062470 -0.587196 -0.724353 -1.550162 1.510451 0.304944 > 0.024780 -0.202001 0.582708 1.729426 1.843437 0.907140 > 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 > 0.000000 0.000000 0.000000 0.000000 0.000000 > Parameter 10 to routine DGEMM was incorrect > Mac OS BLAS parameter error in DGEMM , parameter #0, (unavailable), is 0 > > and R immediately dies. I know the arguments are being passed into > the C code and everything up to my F77_CALL is functioning based on > the printed output. The problem is definitely something to do with my > F77_CALL(dgemm) line. My understanding is that parameter 10 should be > the leading dimension of the matrix B, which in this case should be > equal to 2, the number of rows in that matrix, which is what I am > doing. I have also considered that parameter numbering starts at 0, > in which case the incorrect parameter is &zero, but again that seems > correct to me. All of my reading and research suggests I am doing > everything correctly, so I am somewhat stumped. Perhaps I am missing > something simple or obvious, as I have never done this before and am > proceeding with only google and the R docs as my guide. I am > wondering if anybody can see what I'm doing wrong here, or perhaps > something I could do to try to fix it. Any assistance would be > greatly appreciated. > > Best Regards, > > Jason Rudy > Graduate Student > Bioinformatics and Medical Informatics Program > San Diego State University > > ______________________________________________ > R-devel at r-project.org mailing list > https://stat.ethz.ch/mailman/listinfo/r-devel >-- Brian D. Ripley, ripley at stats.ox.ac.uk Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/ University of Oxford, Tel: +44 1865 272861 (self) 1 South Parks Road, +44 1865 272866 (PA) Oxford OX1 3TG, UK Fax: +44 1865 272595
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