[Rd] Using matprod from array.c

[Heather Turner]

Hi,

I was wondering if I could use the matprod function from array.c in my own C routine.  I tried the following as a test

/* my_matprod.c */

# include <Rinternals.h> /* for REAL, SEXP etc */
# include <R_ext/Applic.h> /* array.c says need for dgemm */

/* following copied from array.c */
static void matprod(double *x, int nrx, int ncx,
		    double *y, int nry, int ncy, double *z)
{
    char *transa = "N", *transb = "N";
    int i,  j, k;
    double one = 1.0, zero = 0.0, sum;
    Rboolean have_na = FALSE;

    if (nrx > 0 && ncx > 0 && nry > 0 && ncy > 0) {
	/* Don't trust the BLAS to handle NA/NaNs correctly: PR#4582
	 * The test is only O(n) here
	 */
	for (i = 0; i < nrx*ncx; i++)
	    if (ISNAN(x[i])) {have_na = TRUE; break;}
	if (!have_na)
	    for (i = 0; i < nry*ncy; i++)
		if (ISNAN(y[i])) {have_na = TRUE; break;}
	if (have_na) {
	    for (i = 0; i < nrx; i++)
		for (k = 0; k < ncy; k++) {
		    sum = 0.0;
		    for (j = 0; j < ncx; j++)
			sum += x[i + j * nrx] * y[j + k * nry];
		    z[i + k * nrx] = sum;
		}
	} else
	    F77_CALL(dgemm)(transa, transb, &nrx, &ncy, &ncx, &one,
			    x, &nrx, y, &nry, &zero, z, &nrx);
    } else /* zero-extent operations should return zeroes */
	for(i = 0; i < nrx*ncy; i++) z[i] = 0;
}

/* test function: matrix multiplication of nr by nc matrix with nc-vector */
SEXP my_matprod(SEXP M, SEXP v, SEXP nr, SEXP nc) {
  R_len_t  nrm = INTEGER(nr)[0], ncm = INTEGER(nc)[0];
  SEXP ans;

  PROTECT(ans = allocMatrix(REALSXP, nrm, 1));
  matprod(REAL(M), nrm, ncm,
	  REAL(v), ncm, 1, REAL(ans));
  UNPROTECT(1);
  return(ans);
}

When I try to make the DLL I get the following
D:\C_routines>RCMD SHLIB my_matprod.c
making my_matprod.d from my_matprod.c
gcc   -IC:/R/tex/R-2.2.0/include -Wall -O2   -c my_matprod.c -o my_matprod.o
ar cr my_matprod.a my_matprod.o
ranlib my_matprod.a
gcc  --shared -s  -o my_matprod.dll my_matprod.def my_matprod.a  -LC:/R/tex/R-2.
2.0/src/gnuwin32   -lg2c -lR
my_matprod.a(my_matprod.o)(.text+0x19a):my_matprod.c: undefined reference to `dg
emm_'
collect2: ld returned 1 exit status
make: *** [my_matprod.dll] Error 1

I'm using Windows XP and I'm using the MinGW gcc. It works fine if I comment out the Fortran call so that the method for have_na is always used. 

Do I need to include another header file to use dgemm? Is there a better way to use matprod than just to copy the code?

Any help appreciated,

Heather


Dr H Turner
Research Assistant
Dept. of Statistics
The University of Warwick
Coventry
CV4 7AL

Tel: 024 76575870
Url: www.warwick.ac.uk/go/heatherturner