[Help-gsl] Seg faults and altered parameters using vector views

[Michael Braun]

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