From ripley at stats.ox.ac.uk  Mon Feb 21 08:42:00 2011
From: ripley at stats.ox.ac.uk (Prof Brian Ripley)
Date: Mon, 21 Feb 2011 07:42:00 +0000 (GMT)
Subject: [Rd] Problem using F77_CALL(dgemm) in a package
In-Reply-To: <D11E54DC-6FED-4AE3-A5DB-7EDBDBA28DD4@r-project.org>
References: <AANLkTimAHcu_=_rMR9bFVYK86bmWhs8Pimgv3L-LiHsV@mail.gmail.com>
	<alpine.LFD.2.02.1102201540060.16292@gannet.stats.ox.ac.uk>
	<AANLkTi=CQhxm5RBhymbCBLz8JAYs1cbEy=xC9LR_n8js@mail.gmail.com>
	<D11E54DC-6FED-4AE3-A5DB-7EDBDBA28DD4@r-project.org>
Message-ID: <alpine.LFD.2.02.1102210738150.31556@gannet.stats.ox.ac.uk>

Simon,

FWIW, %*% does a bit more, in particular does not call dgemm with NAs 
present, as BLAS are often 'optimized' to give the wrong answer in 
that case (e.g. by assuming 0*x is always 0, even though x can be 
Inf, or my not distinguishing NaNs, whereas R uses one for NA).

Brian

On Sun, 20 Feb 2011, Simon Urbanek wrote:

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



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