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@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
