[Top][All Lists]
[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[Help-gsl] Seg faults and altered parameters using vector views
From: |
Michael Braun |
Subject: |
[Help-gsl] Seg faults and altered parameters using vector views |
Date: |
Thu, 4 Jan 2007 17:59:44 -0500 |
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
address@hidden
(617) 253-3436
- [Help-gsl] Seg faults and altered parameters using vector views,
Michael Braun <=