I am having some trouble working with vector views, and I am hoping someone can point me in the right direction (I am a relatively new C programmer, so I am still learning some new tricks).
Below is a function that computes the log density of a multivariate normal distribution using gsl vectors and matrices. The eventual result is correct (i.e., consistent with other sources). The problem is that the matrix sigma is changed in the calling function (the top left corner element changes to something else). I find this odd, since I never change any part of sigma in this function. In addition, I draw your attention to the two fprintf statements in the lines marked A and B. If both are commented out, there is no problem other than what I just described. If the fprintf statement is inserted at line A *or before*, again there is no additional problem. However, if I call the fprintf function at line B or later, I get a segmentation fault with the message *** caught segfault *** address 0x4, cause 'memory not mapped' (the 4 in this statement is the value for k; if I change k, the address changes as well). So it looks like I am doing something wrong with pointers and memory allocation, but I don't know what. I still get the right "answer" in res, but I'd like to understand this other behavior. Here is the code.. I am calling this code from R, so the following fucntion converts the R objects to GSL vectors and matrices, calls the logpdf function, and returns the value to R: SEXP R_mvnorm_logpdf (SEXP xx, SEXP mux, SEXP sigmax, SEXP kx) { // These lines copy the R objects (since they are not copied in the R call) SEXP sigmaCopy, xCopy, muCopy, res; PROTECT(sigmaCopy = duplicate(sigmax)); PROTECT(xCopy = duplicate(xx)); PROTECT(muCopy = duplicate(mux)); // Pointers to the "data" in the SEXP objects double * xAr = REAL(xCopy); double * muAr = REAL(muCopy); double * sigmaAr = REAL(sigmaCopy); // GSL views from the arrays in the R SEXP objects 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; // call logpdf function double logans = gsl_MB_mvnorm_logpdf(x, mu, sigma, k); // return result to R (this works fine) PROTECT(res=allocVector(REALSXP,1)); REAL(res)[0] = logans; UNPROTECT(4); return(res); } This is the function that computes the log density (and this is where the seg fault takes place) double gsl_MB_mvnorm_logpdf(gsl_vector * beta, gsl_vector * betaMean, gsl_matrix * sigma, int k) { // computes density of multivariate normal vector at point 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); 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) A: // gsl_matrix_fprintf(stdout,sigma,"%f"); gsl_vector * x = gsl_vector_alloc(k); B: // 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 I did notice a thread in the help-gsl archives that discusses how to use views, but none of the proposed solutions worked for me. So I am flummoxed, and I would greatly appreciate any help you can provide. Best wishes, Michael Braun ------------------------------------------ Michael Braun Assistant Professor of Marketing MIT Sloan School of Management 38 Memorial Drive, E56-329 Cambridge, MA 02139 [EMAIL PROTECTED] (617) 253-3436 _______________________________________________ Help-gsl mailing list Help-gsl@gnu.org http://lists.gnu.org/mailman/listinfo/help-gsl