[R] Problem with .C()

[Rolf Turner]

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}}