[Rd] Parameter changes and segfault when calling C code through .Call

[Michael Braun]

I am experiencing some odd behavior with the .Call interface, and I am hoping someone out there can help me with it.  Below is a simple example (in that there are R packages that do exactly what I want), but this code illustrates the problem.  Thanks in advance for any help you can provide.

Suppose I want to compute the log density of a multivariate normal distribution using C code and the gsl library.  My R program is:

	dyn.load("mvnorm-logpdf.so")

	x<-c(0,0,0,0,0,0)
	mu<-c(0,0,0,0,0,0)
	sig<-diag(6)
	print(sig)
	w<-.Call("R_mvnorm_logpdf",as.double(x),as.double(mu),sig, as.integer(6))
	print(sig) # sig has changed after .Call

This code takes the SEXP's that were passed from R and converts them to gsl objects, and then calls the logpdf function:

# include <R.h>
# include <Rinternals.h>
# include <Rmath.h>
# include <gsl/gsl_matrix.h>
# include <gsl/gsl_vector.h>
# include <gsl/gsl_blas.h>
# include <gsl/gsl_linalg.h>
# include <gsl/gsl_math.h>


	SEXP R_mvnorm_logpdf (SEXP xx, SEXP mux, SEXP sigmax, SEXP kx) {
	
		int k = INTEGER(kx)[0];
		double * xAr = REAL(xx);
		double * muAr = REAL(mux);
		double * sigmaAr = REAL(sigmax);
		SEXP res;
	
		gsl_vector_view xView = gsl_vector_view_array(xAr,k);
		gsl_vector_view muView = gsl_vector_view_array(muAr,k);
		gsl_matrix_view sigmaView = gsl_matrix_view_array(sigmaAr,k,k);

		gsl_vector * x = &xView.vector;
		gsl_vector * mu = &muView.vector;
		gsl_matrix * sigma = &sigmaView.matrix;

1:		double logans = gsl_MB_mvnorm_logpdf(x, mu, sigma, k); // <-call logpdf here

		PROTECT(res=allocVector(REALSXP,1));
		REAL(res)[0] = logans;
		UNPROTECT(1);
		return(res);
	}

The logpdf function is here

	double gsl_MB_mvnorm_logpdf(gsl_vector * beta, gsl_vector * betaMean, gsl_matrix * sigma, int k) {
		// computes density of multivariate normal vector at vector beta, with mean betaMean and cov sigma

		double logdetSigma = 0;
		double res;
		double * kern;
		int i, err;

		// pointer to Cholesky decomp of sigma
		gsl_matrix * sigmaChol = gsl_matrix_alloc(k, k); // define matrix that will store Chol decomp
		gsl_matrix_memcpy(sigmaChol, sigma);
		gsl_linalg_cholesky_decomp(sigmaChol);

		// compute logdet of sigma by 2*sum of log of diagomal elements of chol decomp
		for (i=0; i<k; i++) {
			logdetSigma = logdetSigma + log(gsl_matrix_get(sigmaChol,i,i));
		}
		logdetSigma = 2*logdetSigma;
	
		// compute (beta-mean)' sigma^(-1) (beta-mean)
		gsl_vector * x = gsl_vector_alloc(k);

2:  // gsl_matrix_fprintf(stdout,sigma,"%f");

		gsl_vector_memcpy(x, beta);
		gsl_vector_sub(x, betaMean); // beta - betaMean
	
		gsl_vector * y = gsl_vector_alloc(k);
		gsl_vector_memcpy(y,x);
		gsl_blas_dtrsv(CblasLower, CblasNoTrans, CblasNonUnit, sigmaChol, y); // y = inv(chol)*x from BLAS
		gsl_blas_ddot(y,y,kern); // kern = y'y

		// compute log density
		res = -k*M_LN_SQRT_2PI - 0.5*(logdetSigma + *kern);
 
		// release space
		gsl_matrix_free(sigmaChol);
		gsl_vector_free(x);
		gsl_vector_free(y);

		return(res);
	} // end gsl_mvnorm_pdf

The problem is that after I make the .Call in R, the value of the sig matrix changes (the 1 in the upper left corner changes from a 1 to a 0).  Since I don't make changes to the sigma object directly, I don't know how that could happen.  The following output is the sig matrix, before and after the .Call.

	> source("pdfTest.R")
	     [,1] [,2] [,3] [,4] [,5] [,6]
	[1,]    1    0    0    0    0    0
	[2,]    0    1    0    0    0    0
	[3,]    0    0    1    0    0    0
	[4,]    0    0    0    1    0    0
	[5,]    0    0    0    0    1    0
	[6,]    0    0    0    0    0    1
	     [,1] [,2] [,3] [,4] [,5] [,6]
	[1,]    0    0    0    0    0    0
	[2,]    0    1    0    0    0    0
	[3,]    0    0    1    0    0    0
	[4,]    0    0    0    1    0    0
	[5,]    0    0    0    0    1    0
	[6,]    0    0    0    0    0    1

Through some debugging, I do know that it occurs somewhere in the logpdf function (line 1: in the first func).  I should also note that this function does compute the density correctly.

Here's some other, potentially relevant information.  Note the print statement in line 2 of the logpdf function.  If that print were "uncommented" and placed somewhere *earlier* in the function, the elements of the sigma matrix are as they should be, and the program runs with no problem (except for the change in sig as described above).  However, if the print statement is at line 2: or later, the correct matrix elements are printed, but then the .Call crashes with the following message:

	 *** caught segfault ***
	address 0x6, cause 'memory not mapped'

(If I change the size of the matrix to diag(k), the 0x6 becomes 0xk). So, for some reason, k is interpreted as a memory address.  I double-checked for places where I may have confused pointers and values, but I can't see where that would have happened.  Also, after searching through the R-devel archives, I noticed that others have had other odd 'memory not mapped' errors in totally different scenarios.

So I am flummoxed, and don't really know where to go from here.

Best wishes,

MB



------------------------------------------
Michael Braun
Assistant Professor of Marketing
MIT Sloan School of Management
38 Memorial Drive, E56-329
Cambridge, MA 02139
braunm at mit.edu
(617) 253-3436