From jcrudy at gmail.com  Sun Feb 20 23:23:26 2011
From: jcrudy at gmail.com (Jason Rudy)
Date: Sun, 20 Feb 2011 14:23:26 -0800
Subject: [Rd] Problem using F77_CALL(dgemm) in a package
In-Reply-To: <alpine.LFD.2.02.1102201540060.16292@gannet.stats.ox.ac.uk>
References: <AANLkTimAHcu_=_rMR9bFVYK86bmWhs8Pimgv3L-LiHsV@mail.gmail.com>
	<alpine.LFD.2.02.1102201540060.16292@gannet.stats.ox.ac.uk>
Message-ID: <AANLkTi=CQhxm5RBhymbCBLz8JAYs1cbEy=xC9LR_n8js@mail.gmail.com>

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
>



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