Dear all,
I need to calculate tr(xyxy) fast for matrices x and y. Inspired by the
R-source code, I've created the following functions (I am *new* to writing
external C-functions, so feel free to laugh at my code - or, perhaps, suggest
changes):
#include <Rinternals.h>
#include <R_ext/Applic.h> /* for dgemm */
static void matprod(double *x, int nrx, int ncx,
double *y, int nry, int ncy, double *z)
{
char *transa = "N", *transb = "N";
double one = 1.0, zero = 0.0;
F77_CALL(dgemm)(transa, transb, &nrx, &ncy, &ncx, &one,
x, &nrx, y, &nry, &zero, z, &nrx);
}
SEXP trProd2(SEXP x, SEXP y)
{
int nrx, ncx, nry, ncy, mode, i;
SEXP xdims, ydims, ans, ans2, tr;
xdims = getAttrib(x, R_DimSymbol);
ydims = getAttrib(y, R_DimSymbol);
mode = REALSXP;
nrx = INTEGER(xdims)[0];
ncx = INTEGER(xdims)[1];
nry = INTEGER(ydims)[0];
ncy = INTEGER(ydims)[1];
PROTECT(ans = allocMatrix(mode, nrx, ncy));
PROTECT(ans2 = allocMatrix(mode, nrx, ncy));
PROTECT(tr = allocVector(mode, 1));
matprod(REAL(x), nrx, ncx, REAL(y), nry, ncy, REAL(ans));
matprod(REAL(ans), nrx, ncy, REAL(ans), nrx, ncy, REAL(ans2));
for (i=0; i< nrx; i++){
REAL(tr)[0] = REAL(tr)[0] + REAL(ans2)[i*(nrx+1)];
}
UNPROTECT(3);
return(tr);
}
In R I do:
trProd2 <- function(x, y) {
.Call("trProd2", x, y)
}
x <- y <- matrix(1:4,ncol=2)
storage.mode(x) <- storage.mode(y) <- "double"
for (i in 1:10)
print(trProd2(x, y))
[1] 835
[1] 833
[1] 833
[1] 862
[1] 834
[1] 835
[1] 834
[1] 836
[1] 862
[1] 833
The correct answer is 833. Can anyone give me a hint to what is wrong? I am on
a windows xp platform.
Thanks in advance
Søren
______________________________________________
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
