I've been trying to get my head around using matrices in calls to .C().
As an exercise I wrote some code to calculate the product of two
matrices.
(Well, it makes it easy to check if one is getting the right answer!)
After obtaining some advice from a Certain Very Wise Person at Oxford,
(to find out how to deal with array indexing in C functions called from
elsewhere) I wrote the following C code:
# define MAT(X,I,J,M) (X[(I)+(M)*(J)])
void matmult(x,y,mx,nx,my,ny,z)
double *x, *y, *z;
int *mx, *nx, *my, *ny;
{
int vmx, vnx, vmy, vny, i, j, k;
vmx = *mx;
vnx = *nx;
vmy = *my;
vny = *ny;
for(i=0; i<vmx; i++) {
for(j=0; j<vny; j++) {
MAT(z,i,j,vmx) = 0.;
for(k=0; k<vny; k++) {
MAT(z,i,j,vmx) = MAT(z,i,j,vmx) + MAT
(x,i,k,vmx)* MAT(y,k,j,vmy);
}
}
}
return;
}
I wrote R code to call this function as follows:
matmult <- function(x,y) {
if(ncol(x) != nrow(y)) stop("Non-conformable dimensions.")
if(!is.loaded("matmult")) dyn.load("matmult.so")
z <- .C(
"matmult",
x=as.double(x),
y=as.double(y),
mx=as.integer(nrow(x)),
nx=as.integer(ncol(x)),
my=as.integer(nrow(y)),
ny=as.integer(ncol(y)),
z=double(nrow(x)*ncol(y))
)$z
matrix(z,nrow=nrow(x),ncol=ncol(y))
}
I did the R CMD SHLIB thing, started R, sourced in the code for the
matmult R function,
created matrices
a <- matrix(1:9,3,3)
b <- matrix(1:15,3,5)
and did
matmult(a,b)
I got
R See > matmult(a,b)
[,1] [,2] [,3] [,4] [,5]
[1,] 2.553124e+302 4.084999e+302 5.616873e+302 7.148748e+302 183.8281
[2,] 3.600000e+01 8.100000e+01 1.260000e+02 1.710000e+02 216.0000
[3,] NaN NaN NaN NaN NaN
Persisting (perversely) I repeated the calculation (changing *nothing*)
several times:
R See > matmult(a,b)
[,1] [,2] [,3] [,4] [,5]
[1,] 3.000000e+01 6.600000e+01 1.020000e+02 1.380000e+02 NaN
[2,] 3.600000e+01 8.100000e+01 1.260000e+02 1.710000e+02 NaN
[3,] 9.474724e+64 1.658077e+65 2.368681e+65 3.079285e+65 NaN
R See > matmult(a,b)
[,1] [,2] [,3] [,4] [,5]
[1,] 1.846886e+20 2.955017e+20 4.063149e+20 5.17128e+20 174
[2,] 3.600000e+01 8.100000e+01 1.260000e+02 1.71000e+02 216
[3,] NaN NaN NaN NaN NaN
R See > matmult(a,b)
[,1] [,2] [,3] [,4] [,5]
[1,] -1.355635e+69 -2.372361e+69 -3.389086e+69 -4.405812e+69 174
[2,] -3.564267e+292 -6.237466e+292 -8.910666e+292 -1.158387e+293 216
[3,] NaN NaN NaN NaN NaN
R See > matmult(a,b)
[,1] [,2] [,3] [,4] [,5]
[1,] -2.682205e+283 -4.693860e+283 -6.705514e+283 -8.717168e+283 174
[2,] 2.033693e+164 3.558964e+164 5.084234e+164 6.609504e+164 216
[3,] NaN NaN NaN NaN NaN
R See > matmult(a,b)
[,1] [,2] [,3] [,4] [,5]
[1,] 3.000000e+01 6.600000e+01 1.020000e+02 1.380000e+02 174
[2,] 6.055777e+111 9.689243e+111 1.332271e+112 1.695618e+112 216
[3,] 4.200000e+01 9.600000e+01 1.500000e+02 2.040000e+02 258
R See > matmult(a,b)
[,1] [,2] [,3] [,4] [,5]
[1,] -1.355635e+69 -2.372361e+69 -3.389086e+69 -4.405812e+69 174
[2,] -3.564267e+292 -6.237466e+292 -8.910666e+292 -1.158387e+293 216
[3,] 4.200000e+01 9.600000e+01 1.500000e+02 2.040000e+02 258
R See > matmult(a,b)
[,1] [,2] [,3] [,4] [,5]
[1,] -6.758072e+304 -1.182663e+305 -1.689518e+305 -2.196373e+305 174
[2,] 6.629826e+179 1.160220e+180 1.657457e+180 2.154694e+180 216
[3,] -5.288177e+276 -9.254310e+276 -1.322044e+277 -1.718658e+277 258
R See > matmult(a,b)
[,1] [,2] [,3] [,4] [,5]
[1,] 30 66 102 138 174
[2,] 36 81 126 171 216
[3,] 42 96 150 204 258
this last result being the right answer.
Several more repetitions gave the right answer, over and over again,
and then ....
R See > matmult(a,b)
[,1] [,2] [,3] [,4] [,5]
[1,] 3.000000e+01 6.600000e+01 1.020000e+02 1.380000e+02 174
[2,] 2.633743e+171 4.213988e+171 5.794234e+171 7.374479e+171 216
[3,] 4.200000e+01 9.600000e+01 1.500000e+02 2.040000e+02 258
back to rubbish again! (Then another piece of rubbish, then another
right answer ....)
I cannot for the life of me figure out what's going on. Can anyone
out there spot the loony?
I must be doing something dumb, *somewhere*, but I can't see it. The
fact that
I'm getting different answers on different occasions suggests a
memory leak somewhere.
Or a failure to initialize something ... but where?
With eternal gratitude for any enlightenment,
cheers,
Rolf Turner
######################################################################
Attention:\ This e-mail message is privileged and confid...{{dropped:9}}
______________________________________________
[email protected] mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
