Problem using F77_CALL(dgemm) in a package

Jason,

FWIW the direct interface (.Call) is more efficient and makes passing things from R simpler:

C_matrix_multiply = function(A,B) .Call("R_matrix_multiply", A, B)

The drawback is a bit more legwork on the C side, but it also gives you more flexibility:

SEXP R_matrix_multiply(SEXP A, SEXP B) {
? ? ? ?double one = 1.0;
? ? ? ?double zero = 0.0;
? ? ? ?int *dimA = INTEGER(getAttrib(A, R_DimSymbol));
? ? ? ?int *dimB = INTEGER(getAttrib(B, R_DimSymbol));
? ? ? ?SEXP sDimC = PROTECT(allocVector(INTSXP, 2));
? ? ? ?int *dimC = INTEGER(sDimC);
? ? ? ?SEXP C = PROTECT(allocVector(REALSXP, dimA[0] * dimB[1]));
? ? ? ?if (dimA[1] != dimB[0]) error("incompatible matrices!");
? ? ? ?dimC[0] = dimA[0];
? ? ? ?dimC[1] = dimB[1];
? ? ? ?setAttrib(C, R_DimSymbol, sDimC);
? ? ? ?A = PROTECT(coerceVector(A, REALSXP));
? ? ? ?B = PROTECT(coerceVector(B, REALSXP));
? ? ? ?F77_CALL(dgemm)("N","N",dimA,dimB+1,dimA+1,&one,REAL(A),dimA,REAL(B),dimA+1,&zero,REAL(C),dimA);
? ? ? ?UNPROTECT(4);
? ? ? ?return C;
}

For comparison:

A=matrix(rnorm(1e5),500)
B=matrix(rnorm(1e5),,500)

.Call:

system.time(for (i in 1:10) C_matrix_multiply(A,B))

? user ?system elapsed
?0.656 ? 0.008 ? 0.686

.C:

system.time(for (i in 1:10) CC_matrix_multiply(A,B))

? user ?system elapsed
?0.886 ? 0.044 ? 0.943


in fact .Call is even a tiny bit faster than %*%:

system.time(for (i in 1:10) A %*% B)

? user ?system elapsed
?0.658 ? 0.004 ? 0.665

(it's not just a measurement error - it's consistent for more replications etc. - but it's really negligible - possibly just due to dispatch of %*%)

Cheers,
Simon


On Feb 20, 2011, at 5:23 PM, Jason Rudy wrote:

It was indeed a simple problem! ?I took a look at that array.c as you
suggested and that cleared it right up. ?So, the correct C code is:

#include <R.h>
#include <R_ext/Utils.h>
#include <R_ext/Lapack.h>
#include <R_ext/BLAS.h>

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,p,n,&one,A,m,B,n,&zero,C,m);
}

The only difference being that I had the 4th and 5th arguments (n and
p) mixed up. ?There was also a problem in my R code after the
multiplication took place. ?For the record, the correct R code is:

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[[6]],nrow(A),ncol(B)))
}

Thanks for the help. ?Now that I have a functioning example I am well
on my way to completing this project.

-Jason

On Sun, Feb 20, 2011 at 7:42 AM, Prof Brian Ripley
<ripley at stats.ox.ac.uk> wrote:

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

--
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