[Rcpp-devel] The problem i encountered about RcppEigen double -(k/2)*std::log(2*M_PI)+rooti.diagonal().log().sum(); can not run .

[李晓辉]

These days i rewrote the  dmvnorm() function from the mvtnorm package with
RcppEigen which is publised in
http://gallery.rcpp.org/articles/dmvnorm_arma/.
The code as follows:

> #include <RcppEigen.h>
> using namespace Rcpp;
> // [[Rcpp::depends(RcppEigen)]]
> //[[Rcpp::export]]
> using Eigen::Map;
> using Eigen::MatrixXd;
> using Eigen::VectorXd;
> using Eigen::LLT;
> VectorXd dmvnrm_Eigen(Map<MatrixXd> X,Map<VectorXd> mu,Map<MatrixXd>
> Sigma){
>   const int k = X.cols();
>    int n = X.rows();
>    MatrixXd D(k,k);
>    LLT<MatrixXd>lltOfsigma(Sigma);
>    MatrixXd L = lltOfsigma.matrixL();
>    MatrixXd rooti = L.inverse()*Eigen::MatrixXd::Identity(k,k);
>    MatrixXd MU = mu.replicate(1,n);
>    MatrixXd Q = rooti.transpose()*(X.transpose()-MU);
>    MatrixXd QUAD = Q.cwiseProduct(Q);
>    VectorXd quads = QUAD.colwise().sum();
>    VectorXd FD(n);
>    double Con = -(k/2)*std::log(2*M_PI)+rooti.diagonal().log().sum();
>    FD.setConstant(Con);
>    return (FD - 0.5*quads).exp();
> }

/*** R
> sigma <- bayesm::rwishart(10,diag(8))$IW
> means <- rnorm(8)
> X     <- mvtnorm::rmvnorm(90, means, sigma)
> dmvnrm_Eigen(X,means,sigma)
> */


When i run above code with RStudio, it often got an error message *:R
Session Aborted*
 it seems that* double Con =
-(k/2)*std::log(2*M_PI)+rooti.diagonal().log().sum(); *got a
problem.Actually,if i rewrite these code as follow:

> double Con=0.;
>    for (int i=0;i<k;i++){
>    Con+=log(rooti(i,i));
>    }
> Con+=-(k/2)*std::log(2*M_PI)


it work well.
I do not know what is wrong with
*double Con = -(k/2)*std::log(2*M_PI)+rooti.diagonal().log().sum();*
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.r-forge.r-project.org/pipermail/rcpp-devel/attachments/20141029/61339fd0/attachment.html>